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1.  Introduction 


■ 


Consider  the  NxN  system  of  linear  equations 

' (1)  M X - b, 

where  the  coefficient  matrix  M is  large,  sparse,  and  nonsymmetric . 

Assume  that  M can  be  factored  in  the  form 

M - L D U, 

where  L is  a lowr  triangular  matrix,  D is  a diagonal  matrix,  and  U is  a 
unit  upper  triangular  matrix.  Such  systems  arise  frequently  in 
scientific  computation,  e.g.,  in  finite  difference  and  finite  element 
approximations  to  non-self-adjoint  elliptic  boundary  value  problems,  inv 
this  report»-«e  present^ a package  of  efficient,  reliable, 
well-documented,  and  portable  FORTRAN  subroutines  for  solving  these 
systems.,  See  [3]  for  a corresponding  package  for  symmetric  problems. 

Direct  methods  for  solving  (I)  are  generally  variations  of  Gaussian 
elimination.  We  form  the  LDU  decomposition  of  A,  and  successively  solve 
the  triangular  systems 

(2)  L y - b,  D z - y,  U x - z. 

When  M is  large  (N  » 1),  (dense)  Gaussian  elimination  is  prohibitively 

3 

expensive  in  terms  of  both  the  work  ('  2/3  N multiplies)  and  storage 
2 

(N  words)  required.  But,  since  M is  sparse,  most  entries  of  N,  L,  and 
U are  zero  and  there  are  significant  advantages  to  factoring  M without 


storing  or  operating  on  these  xeroes.  Recentlyt  • nissber  of 
Implementations  of  sparse  Gaussian  elimination  have  appeared  based  on 
this  Idea,  c£.,  (2,  6,  7,  8]. 

In  section  2,  we  describe  the  scheme  used  for  storing  sparse 
matrices,  while.  In  section  3,  we  give  an  overview  of  the  package  from 
the  point  of  view  of  the  user;  for  further  details  of  the  algorithms 
employed,  see  [4,  5].  In  section  4,  we  Illustrate  the  performance  of 
the  package  on  a typical  model  problem.  Listings  of  the  three  sets  of 
subroutines  for  factoring  and  solving  the  class  of  sparse  nonsymmetrlc 
systems  under  consideration  appear  In  Appendices  1,  2,  and  3.  These 
three  sets  of  subroutines  have  different  storage  schemes  and  basically 
crade~off  run-time  efficiency  for  storage.  Appendix  4 contains  a test 
driver  which  sets  up  a problem  and  calls  all  three  sets  of  subroutines 
for  solution.  A sample  output  appears  as  Appendix  5. 
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2.  Sparse  Matrix  Storage  Schemes 

Since  the  coefficient  matrix  M and  the  triangular  factors  L and  U 
are  large  and  sparse.  It  Is  Inefficient  to  store  them  as  dense  matrices. 
The  package  has  two  schemes  for  storing  sparse  matrices,  called  the 
"uncompressed  storage  scheme”  and  the  "compressed  storage  scheme."  The 
Input  matrix  M Is  always  stored  using  the  first  of  these,  while  the 
triangular  factors  L and  U may  be  stored  using  either  one,  depending  on 
which  subroutines  are  used.  The  subroutine  NDRV  uses  the  "uncompressed 
storage  scheme"  for  L and  U while  the  subroutines  TDRV  and  CDRV  use  the 
"compressed  storage  scheme." 

The  uncompressed  storage  scheme  has  been  used  previously  In  various 
forms,  cf.  (1,  6].  To  use  It  to  store  the  Input  matrix  H requires 
three  one-dlmenslonal  arrays:  lA,  JA,  and  A.  The  nonzero  en  les  of  M 
are  stored  row-by-row  in  the  REAL  array  A.  To  identify  the  individual 
nonzero  entries  In  a row,  we  need  to  know  in  which  column  each  entry 
lies.  The  INTEGER  array  JA  contains  the  column  Indices  which  correspond 
to  the  nonzero  entries  of  M,  l.e..  If  A(K)  ■ M(1,J),  then  JA(K)  > J.  In 
addition,  we  need  to  know  where  each  row  starts  and  how  long  It  Is.  The 
INTEGER  array  lA  contains  the  Index  positions  In  JA  and  A where  the  rows 
of  M begin,  l.e..  If  M(I,J)  is  the  first  (leftmost)  entry  of  the  I-th 
row  and  A(K)  - M(I,J),  then  IA(I)  - K.  Moreover,  IA(N+1)  Is  defined  as 
the  Index  In  JA  and  A of  the  first  location  following  the  last  element 
In  the  last  row.  Thus,  the  nu,aber  of  entries  in  the  I-th  row  is  given 
by  IA(I+1)  - IA(I),  the  nonzero  entries  of  the  I-th  row  are  stored 


consecutively  in 


A(1A(I)),  A(IA(l)+l) A(lA(I-*-l)-l), 

and  the  corresponding  colmn  indices  are  stored  consecutively  in 
JA(IA(I)),  JA(IA(1)+1) JAdAd+D-l). 

For  example,  the  SxS  matrix 

1.  0.  2.  0.  0. 

0.  3.  0.  0.  0. 

M - 0.  4.  3.  6.  0. 

0.  0.  0.  7.  0. 

0.  0.  0.  8.  9. 

is  stored  as 

11  2 3 4 5 6 7 8 9 

lA  I 1 3 4 7 8 10 

JAI1  32  23  44  45 

A I 1.  2.  3.  4,  5.  6.  7.  8.  9.  . 


The  overhead  in  this  storage  scheme  is  the  storage  required  for  the 
INTEGER  arrays  lA  and  JA.  But  since  lA  has  N+1  entries  and  JA  has  one 
entry  for  each  element  of  A,  the  total  overhead  is  approximately  equal 


to  the  number  of  nonzero  entries  in  M 


^ 

- 6 - 

The  triangular  laacrlces  L and  U are  stored  In  basically  the  saae 
fashion  using  the  arrays  IL,  JL,  L and  ID,  JU,  U respectively t except 
that  the  diagonal  entries  are  not  stored  In  these  arrays.  The  diagonal 
entries  of  L and  U are  known  to  be  ones  and  are  not  stored  and  the  array 
D Is  used  to  store  the  reciprocals  of  the  diagonal  entries  of  the 
diagonal  matrix  D. 

In  certain  situations,  where  storage  Is  at  a premium,  it  Is 
essential  to  reduce  storage  overhead,  even  at  the  cost  of  decreased 
runtime  efficiency.  This  can  be  done  by  storing  L and  U with  the  more 
complex  compressed  storage  scheme.  This  scheme  Incurs  more  operational 
overhead  than  the  uncompressed  storage  scheme,  but  In  many  Important 
cases  the  storage  requirement  can  be  substantially  reduced.  For  a 
detailed  description,  see  (4,  5,  91. 

[ 


1 


r 
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3.  A Sparse  Nonaymmetrlc  Matrix  Package 

The  package  consists  of  a test  driver,  three  driver  subroutines, 
and  eight  subroutines  (see  Figure  1).  The  three  drivers  (subroutines 
NDRV,  TDRV,  and  CORV)  are  specific  impleaentatlon  designs  which 
illustrate  the  space-time  tradeoff  mentioned  in  section  2.  The  test 
driver  (subroutine  NSTST)  sets  up  a model  sparse  nonsyametric  system  of 
linear  equations  and  calls  each  of  the  three  driver  subroutines  to  solve 
the  linear  system.  In  the  remainder  of  this  section,  we  describe  each 
of  these  subroutines  in  somewhat  greater  detail.  The  codes  themselves 
are  extensively  documented;  for  further  details  about  the  algorithms 
employed  see  [4,  5]. 


Figure  1:  A scheaaclc  overview  of  the  aperee  aymmatrlc  eatrlx  package 


I NSTST  I 


I NSF  I I NNF  I I NNS  I f TRK  | | NROC  | | NSFC  | | NMFC  | | NNSC  | 
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Our  basic  design  for  the  iapleaencation  of  sparse  eliainacion 
follows  chat  of  Chang  [1],  which  has  proved  to  be  especially  robust. 

The  first  implementation,  NDRV,  is  designed  for  speed.  It  uses  the 
uncompressed  storage  scheme  for  M,  L,  and  U because  of  the  smaller 
operational  overhead  associated  with  it.  We  break  the  computation  up 
into  three  distinct  phases:  symbolic  factorization  (subroutine  NSF), 
numeric  factorization  and  the  solution  for  one  right-hand  side 
(subroutine  NNF) , and  forward  and  back  solution  for  additional 
right-hand  sides  (subroutine  NNS).  The  subroutine  NSF  computes  the  zero 
structures  of  L and  U from  that  of  M (disregarding  the  nmerical  entries 
in  M) . The  subroutine  NNF  then  uses  the  structural  information 
generated  by  NSF  to  compute  the  numerical  entries  of  L and  U and  to 
solve  for  one  right-hand  side. 

The  main  advantage  of  splitting  up  the  computation  in  this  way  is 
flexibility.  To  solve  a single  system  of  equations,  it  suffices  to  use 
NSF  and  NNF  (PATH>1  in  NDRV).  It  should  be  pointed  out  here  that  a one 
line  modification  of  NNF  can  be  made  to  allow  the  solution  of  a single 
system  without  storing  L:  simply  comment  out  the  line 

L(l)  - - LI, 

as  indicated  in  the  code.  This  change  will  yield  substantial  storage 
savings  without  the  loss  in  efficiency  incurred  by  TDRV.  To  solve 
several  systems  in  which  the  coefficient  matrices  have  the  same  zero 
structure,  it  suffices  to  use  NSF  and  NNF  only  once  each  for  the  first 
system  and  then  to  use  NNF  once  for  each  subsequent  system  (PATH~2  in 
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NDRV).  Finally,  to  solve  several  systems  with  the  same  coefficient 
matrix  but  different  right  hand  sides,  it  suffices  to  use  NSF  and  NNF 
only  once  each  for  the  first  system,  and  then  to  use  NNS  once  for  each 
subsequent  system  (PATH-3  in  NDRV). 

A drawback  to  the  multi-phase  design  of  NDRV  is  that  it  is 
necessary  to  store  the  description  of  the  zero  structures  of  both  L and 
U.  By  giving  up  some  flexibility,  the  second  implementation,  TDRV, 
greatly  reduces  the  storage  requirements.  The  entire  computation  is 
performed  in  a single  phase  (subroutine  TRK)  to  avoid  storing  either  the 
description  or  the  numerical  entries  of  L.  Moreover,  U is  stored  with 
the  compressed  storage  scheme  to  reduce  the  storage  overhead.  This 
subroutine  Incurs  more  operational  overhead  than  NDRV,  and  we  lose  the 
ability  to  efficiently  solve  a sequence  of  related  systems.  However  the 
total  storage  requirements  are  usually  significantly  smaller  (see  Tables 
4-6). 

Finally,  the  third  Implementation,  CDRV,  attempts  to  balance  the 
design  goals  of  speed,  flexibility,  and  storage  economy.  It  splits  the 
computation  as  in  NDRV  to  allow  flexibility  and  efficiency,  but  it  uses 
the  compressed  storage  scheme  as  in  TDRV  to  reduce  storage  overhead. 

The  rows  and  columns  of  the  original  matrix  M can  be  reordered  (e.g.,  to 
reduce  fillln  or  ensure  numerical  stability)  before  calling  CDRV.  If  no 
reordering  is  done,  then  set  R(I)-C(I)“IC(I)  - I for  I-1,...,N.  The 
solution  Z is  returned  in  the  original  order.  If  the  columns  have  been 


reordered  (l.e.,  C(I).NE.I  for  some  I),  then  CDRV  will  call  a subroutine 
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NROC  which  rearranges  each  row  of  JA  and  A.  leaving  the  rows  In  cite 
^ orglnal  order,  but  placing  Che  eleaents  of  each  row  In  Increasing  order 

with  respect  Co  Che  new  ordering.  If  PATH.NE. 1,  then  NROC  Is  assuaed  to 
have  been  called  already. 

To  solve  a single  syccem  of  equations.  It  suffices  to  use  NROC  (If 
Che  columns  of  M have  been  reercered) , NSFC,  and  NNFC  (PATH~1  In  CDRV). 
It  should  be  pointed  out  here  that  a one  line  modification  of  NNFC  can 
be  made  to  allow  the  solution  of  a single  system  without  storing  L: 
simply  comment  out  Che  line 

L(IRL(1))  - - LKl, 

as  Indicated  In  Che  code.  This  change  will  yield  substantial  storage 
savings  without  the  loss  in  efficiency  incurred  in  TDRV.  To  solve 
several  systems  in  which  the  coefficient  matrices  have  the  same  zero 
structure,  it  suffices  to  use  NROC,  NSFC,  and  NNFC  only  once  each  for 
Che  first  system  and  then  to  use  NNFC  once  for  each  subsequent  system 
(PATH-2  In  CORV).  Finally,  to  solve  several  systems  with  the  same 
coefficient  matrix  but  different  right  hand  sides,  it  suffices  to  use 
NR(X:,  NSFC,  and  NNFC  only  once  each  for  the  first  system,  and  Chen  to 
use  NNSC  once  for  each  subsequent  system  (PATH-3  in  CDRV). 

The  test  driver  (program  NSTST)  is  used  to  verify  the  performance 
of  Che  package  on  a particular  computer  system,  and  may  be  used  as  a 
guide  to  understanding  how  to  use  the  package.  It  generates  the 
coefficient  matrix  for  a nonsymmetrlc  five-point  difference  equation  on 


ik. 


- u 


a 3x3  grid  and  chooses  che  right-hand  side  so  that  the  solution  vector  x 
Is  (1 , 2, 3, 4, S, b, 7, b, 9)  (see  Appendix  4).  The  grid  points  are  given  In 
the  natural  row-by-row  ordering.  At  each  stage  the  values  of  all 
relevant  variables  are  printed  out,  and  a sample  output  appears  as 
Appendix  5. 
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4.  Performance 

One  of  the  most  Important  aspects  of  any  package  Is  its  performance 
in  terms  of  both  the  time  and  storage  required  to  solve  a typical 
problem.  In  Tables  1-6,  we  present  the  time  and  storage  required  to 
solve  a nonsymmetric  five-point  difference  equation  on  an  nxn  grid  for 
several  values  of  n.  These  computations  were  performed  in  single 
precision  on  an  IBM  370/158  using  the  FORTRAN  IV  Level  H Extended 
compiler. 

Table  1:  Times  for  S-point  operator  on  a 20x20  mesh 


Code 


NSF(C) 


NNF(C) 


sec/* 


NNS(C) 


Total 


NDRV 

1 

1 

1 

1 

0.213  1 

1 

0.560 



1 1 

1 9.978  1 

1 1 

1 

0.063  1 0.773 

1 

TDRV 

1 

1 

1 

1 

1 

1 

1 1 

1 15.561  1 

1 1 

1 

1 0.873 

I 

1 

1 

1 1 

CDRV 

1 

1 

1 

1 

0.267  1 

1 

0.790 

1 1 

1 14.077  1 

1 1 

1 

0.087  1 1.057 

1 

1 

\ 

\ i 

1 
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Table  2:  Times  for  5-polnt  operator  on  a 30x30  mesh 


r 


Table  4:  Storage  for  S-polnt  operator  on  a 20x20  neah 


Code 


A/JA 


JL 


JU 


Total 


Hults. 


NDRV 


1,920 


3,268 


3,368 


3,368 


3,368 


15,47 


36,118 


TDRV 


1,920 


3,368 


1,889 


7,660 


56,118 


CDRV 


1,920 


3,368 


1,889 


3,368 


1,889 


13,716 


56,118 


Table  5:  Storage  for  5-potnt  operator  on  a 30x30  mesh 


Code 

I 1 i 

! +* 

A/JA  L JL  1 U 

1 - 1 - 1 

JU  1 Total  1 Mults. 

1 1 

NDRV 

1 1 1 

4,380  1 9,456  | 9,456  1 9,456 

* 1 1 1 

1 1 ^ ^ -L 

1 1 

9,456  1 42,327  | 215,708 

1 1 

1 1 

TDRV 

1 1 1 

4,380  1 1 1 9,456 

1 1 1 

1 1 . 1 . . 

1 . 1 - . 

1 1 

4,538  1 19,397  | 215,708 

1 1 

1 1 

CDRV 

1 1 1 

4,380  1 9,456  | 4,538  | 9,456 

1 1 1 

1 1 

4,538  1 35,190  | 215,708 

1 1 

I 


+Total  storage  required  by  driver 
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Table  6:  Storage  for  5-polnC  operator  on  a 40x40  mesh 
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Subroutines  for  Solving  Sparse  Nonsynnetrlc  Systems 
of  Linear  Equations  (Uncompressed  Pointer  Storage) 


C***  Subroutine  NDRV 

C***  Driver  for  subroutines  for  solving  sparse  nonsymmetrlc  systems  of 
C linear  equations  (uncompressed  pointer  storage) 

C 

SUBROUTINE  NDRV 

* (N,  R.C.IC,  lA.JA.A,  B.  Z,  NSP, ISP. RSP. ESP,  PATH,  FLAG) 

C 

C PARAMETERS 

C Class  abbreviations  are  — 

C n - INTEGER  variable 

C f - REAL  variable 

C V - supplies  a VALUE  to  the  driver 

C r - returns  a RESULT  from  the  driver 

C 1 - used  INTERNALly  by  the  driver 

C a - ARRAY 

C 

C Class  I Parameter 

C 1 

C I 

C The  nonzero  entries  of  the  coefficient  matrix  M are  stored 

C row-by-row  in  the  array  A.  To  Identify  the  Individual  nonzero 
C entries  In  each  row,  we  need  to  know  In  which  column  each  entry 

C lies.  The  column  indices  which  correspond  to  the  nonzero  entries 

C of  H are  stored  in  the  array  JA;  l.e.,  if  A(K)  ■ M(I,J),  then 
C JA(K)  ■ J.  In  addition,  we  need  to  know  where  each  row  starts  and 

C how  long  It  Is.  The  Index  positions  in  JA  and  A where  the  rows  of 

C H begin  are  stored  in  the  array  LA;  l.e..  If  M(I,J)  Is  the  first 

C nonzero  entry  (stored)  in  the  I-th  row  and  A(K)  ~ M(I,J),  then 

C 1A(I)  - K.  Moreover,  the  Index  in  JA  and  A of  the  first  location 

C following  the  last  element  In  the  last  row  Is  stored  in  lAfN-fl). 

C Thus,  the  number  of  entries  In  the  I-th  row  is  given  by 
C IA(I+l)  - IA(I),  the  nonzero  entries  of  the  I-th  row  are  stored 
C consecutively  In 

C A(IA(1)),  A(1A(I)+1),  ....  A(IA(1+1)-1), 

C and  the  corresponding  column  Indices  are  stored  consecutively  In 

C JA(IA(1)),  JA(IA(1)+1) JA(IA(I+1)-1). 

C For  example,  the  5 by  3 matrix 

C ( 1.  0.  2.  0.  0.) 

C ( 0.  3,  0.  0.  0.) 

C M - ( 0.  4.  5.  5.  0.) 

C ( 0.  0.  0.  7.  0.) 

C ( 0.  0.  0.  8.  9.) 

C would  be  stored  as 

C 1123456789 


A I 1.  2.  3.  4.  5.  6,  7.  8.  9. 

- number  of  variables/equations. 

- nonzero  entries  of  the  coefficient  matrix  H,  stored 

by  rows. 

Size  - number  of  nonzero  entries  In  M. 

- pointers  to  delimit  the  rows  In  A. 


fva 


f ra 


JA 


Slie  - N+1. 

column  numbers  corresponding  to  the  elamants  of  A. 
Size  • size  of  A. 

right-hand  side  b;  B and  L can  the  same  array. 
Size  • N. 

solution  x;  B and  Z can  be  Che  same  array. 

Size  - N. 


The  rows  and  columns  of  the  original  matrix  H can  be 
reordered  (e.g.,  to  reduce  flllln  or  ensure  niaserlcal  stability) 
before  calling  Che  driver.  If  no  reordering  Is  done. 


R(I)  - C(l)  - 1C(I)  - I 
In  Che  original  order. 


for  1-1, ...,N.  The  solution 


Chen  sec 
Z Is  returned 


H. 


IC 


ordering  of  the  rows  of  M. 

Size  - N. 

ordering  of  the  columns  of 
Size  - N. 

inverse  of  the  ordering  of  Che  columns  of  M; 
1C(C(I))  - 1 for  I-1,.,.,N. 

Size  - N. 


1 .e. , 


The  solution  of  the  system  of  linear  equations  Is  divided  Into 
three  stages  — 

NSF  " The  matrix  M Is  processed  symbolically  to  determine  where 
flllln  will  occur  during  the  numeric  factorization. 

NNF  — The  matrix  M Is  factored  numerically  into  the  product  U>U 
of  a unit  lower  triangular  matrix  L,  a diagonal  matrix  D, 
and  a unit  upper  triangular  matrix  U,  and  the  system 
Mx  - b Is  solved. 

NNS  — The  linear  system  Mx  « b Is  solved  using  the  LDU 
factorization  from  NNF. 

For  several  systems  %diose  coefficient  matrices  have  the  same 
nonzero  structure,  NSF  need  be  done  only  once  (for  Che  first 
system);  then  NNF  Is  done  once  for  each  additional  system.  For 
several  systems  with  the  same  coefficient  matrix,  NSF  and  NNF  need 
be  done  only  once  (for  the  first  system);  then  NNS  Is  done  once 
for  each  additional  right-hand  side. 


PATH  - path  specification;  values  and  their  meanings  are  — 

1 perform  NSF  and  NNF. 

2 perform  NNF  only  (NSF  is  assumed  to  have  been 

done  In  a manner  compatible  with  the  storage 
allocation  used  In  the  driver) . 

3 perform  NNS  only  (NSF  and  NNF  are  assumed  to 

have  been  done  in  a manner  compatible  with  Che 
storage  allocation  used  In  the  driver). 


Various  errors  are  detected  by  the  driver  and  the  Individual 
sub  rout Ines. 


FLAG  - error 
0 

N-«C 

2N4K 

3N+K 

4N+1 

5N4K 

6N-fK 


flag;  values  and  their  meanings  are  — 

No  Errors  Detected 
Null  Row  in  A — Row  - K 
Duplicate  Entry  In  A — Row  - K 
Insufficient  Storage  In  NSF  — Row  - K 
Insufficient  Storage  In  NNF 
Null  Pivot  — Row  - K 
Insufficient  Storage  In  NSF  — Row  - K 


I 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


I 


TN-t-l  Insufficient  Storage  In  NNF 
tM-»K  Zero  Pivot  — Row  • K 
lON-fl  Insufficient  Storsge  in  NORV 

llN-i-t  Illegal  PATH  Specification 


Working  storage  is  needed  for  the  factored  fom  of  the  aatrix 
M plus  various  temporary  vectors.  The  arrays  ISP  and  RSP  should  be 
the  same;  Integer  storage  is  allocated  from  the  beginning  of  ISP 
and  real  storage  from  the  end  of  RSP. 


nv  I NSP 

I 

I 

nvira  1 ISP 


fvira 


RSP 


nr 


ESP 


- declared  dimension  of  ISP  and  RSP;  NSP  generally  must 

be  larger  than  SN-fl  -t-  2K  (where  K - (number  of 
nonzero  entries  in  H) ) . 

- integer  working  storage  divided  up  into  various  arrays 

needed  by  the  subroutines;  ISP  and  RSP  should  be 
the  same  array. 

Size  - NSP. 

- real  working  storage  divided  up  into  various  arrays 

needed  by  the  subroutines;  ISP  and  RSP  should  be 
the  same  array. 

Size  - NSP. 

- if  sufficient  storage  was  available  to  perform  the 

symbolic  factorization  (NSF) . then  ESP  is  set  to  the 
amount  of  excess  storage  provided  (negative  If 
insufficient  storage  was  available  to  perform  the 
numeric  factorization  (NNF)). 


INTEGER  R(l),  C(l),  IC(1),  IA(1),  JA(1),  ISP(l),  ESP, 

* PATH,  FLAG,  Q,  IM,  D,  U,  ROW,  TMP,  UMAX 
REAL  A(l),  B(l),  Z(l),  RSP(l) 


IF  (PATH.LT.l  .OR.  PATH.GT.3)  GO  TO  111 
******  Initialize  and  divide  up  temporary  storage 
IL  - I 

lU  - IL  + N+l 
JL  - lU  + N+l 
FLAG  - 0 


C 

C ******  Call  NSF  if  flag  la  set  ************************************ 

IF  (PATH.GT.l)  GO  TO  2 
IM  - NSP  - N 
Q - IM  - (N+l) 

MAX  - Q - JL 

IF  (MAX.LT.O)  GO  TO  110 
JLMAX  - MAX/2 
JUTMP  - JL  + JLMAX 
JUMAX  - q - JUTMP 
CALL  NSF 

* (N,  R,  IC,  lA,  JA, 

* ISP(IL),  ISP(JL),  JLMAX,  ISP(IU),  ISP(JUTMP),  JUMAX, 

* RSP(q),  RSP(IM),  FLAG) 

IF  (FLAG.NE.O)  GO  TO  100 

Q ******  Move  JU  next  to  JL  ***************************************** 

JLMAX  - lSP(IL+«)-l 
JU  - JL  + JLMAX 
JUMAX  - ISP(IU+N)-1 
IF  (JUMAX.LE.O)  GO  TO  2 
DO  1 J-1, JUMAX 

1 ISP(JU+J-1)  - 1SP(JUTMP+J-1) 


c 

C ******  Call  remaining  subroutines  **************** 

2 JLMAX  - 1SP(IL+N)-1 

JU  - JL  + JLMAX 
JUMAX  - lSP(IU-rtl)-l 
L - JU  + JUMAX 

LMAX  - JLMAX 

D - L + LMAX 
U - D + N 
ROW  - NSP  - N 
IMP  - ROW  - N 
UMAX  - IMP  - U 
ESP  - UMAX  - JUMAX 
C 

IF  (PATH.GT.2)  GO  TO  3 
CALL  NNF 

* (N,  R,  C,  IC,  lA,  JA,  A,  Z,  B, 

* ISP(IL),  ISP(JL),  RSP(L),  LMAX,  RSP(D), 

* ISP(IU),  ISP(JU),  RSP(U).  UMAX, 

* RSP(ROW),  RSP(TMP),  FLAG) 

IF  (FLAG.NE.O)  GO  TO  100 
RETURN 

C 

3 CALL  NNS 

* (N,  R,  C, 

* ISP(IL),  ISP(JL),  RSP(L),  RSP(D), 

* ISP(IU),  ISP(JU),  RSP(U), 

* Z,  B,  RSP(TMP)) 

RETURN 

C 

C **  ERROR:  Error  Detected  in  NSF,  NNF,  or  NNS 
100  RETURN 

C **  ERROR:  InsufElclent  Storage 

110  FLAG  - 10*N  + 1 
RETURN 

G **  ERROR:  Illegal  PATH  Specification 

111  FLAG  - 11*N  + 1 
RETURN 

END 


YALE  SPARSE  MATRIX  PACKAGE  - NONSYMMETRIC  CODES 
SOLVING  THE  SYSTEM  OF  EQUATIONS  Mx  - b 
(UNCOMPRESSED  POINTER  STORAGE) 


I.  SUBROUTINE  NAMES 

Subroutine  nanes  are  of  the  form  Nxx  where  — 

(1)  the  first  letter  Is  N for  nonsynoetrlc  natrlces; 

(2)  the  second  letter  Is  either  S for  symbolic  processing  or  N 

for  numeric  processing; 

(3)  Che  third  letter  Is  either  F for  factorization  or  S for 

solution. 


II.  CALLING  SEQUENCES 

The  coefficient  matrix  can  be  processed  by  an  ordering  routine 
(e.g..  to  reduce  flllln  or  ensure  numerical  stability)  before  using 
the  remaining  subroutines.  If  no  reordering  Is  done,  then  set 
R(I)  " C(I)  - IC(I)  « 1 for  I”1,...,N.  The  calling  sequence  Is  — 


(matrix  ordering)) 

(symbolic  factorization  to  determine  where  flllln  will 
occur  during  numeric  factorization) 

(numeric  factorization  into  product  LDU  of  unit  lower 
triangular  matrix  L,  diagonal  matrix  D,  snd  unit  upper 
triangular  matrix  U,  and  solution  of  linear  system) 
(solution  of  linear  system  for  additional  right-hand 
side  using  LDU  factorization  from  NNF) 


III.  STORAGE  OF  SPARSE  MATRICES 

The  nonzero  entries  of  Che  coefficient  matrix  M are  scored 
row-by-row  In  the  array  A.  To  Identify  the  Individual  nonzero 
entries  in  each  row,  we  need  to  know  in  which  column  each  entry 
lies.  The  column  Indices  which  correspond  to  the  nonzero  entries 
of  M are  stored  In  the  array  JA;  i.e..  If  A(K)  “ M(I,J),  then 
JA(K)  - J.  In  addition,  we  need  to  know  where  each  row  starts  and 
how  long  it  is.  The  index  positions  in  JA  and  A where  Che  rows  of 
H begin  are  stored  in  the  array  lA;  i.e.,  if  H(I,J)  is  the  first 
nonzero  entry  (stored)  in  the  I-th  row  and  A(K)  • M(I,J),  then 
IA(I)  - K.  Moreover,  the  index  in  JA  and  A of  the  first  location 
following  the  last  element  in  the  last  row  Is  stored  In  lA(N-Pl). 
Thus,  the  number  of  entries  In  the  I-th  row  Is  given  by 
IA(I+1)  - IA(I),  the  nonzero  entries  of  the  I-th  row  are  stored 
consecutively  In 

A(IA(1)),  A(IA(I)+1),  ....  A(IA(I+l)-l), 

and  the  corresponding  column  Indices  are  stored  consecutively  in 

JA(1A(I)),  JA(1A(I)+1),  ...,  JA(IA(I+1)-1). 

For  example,  the  3 by  5 matrix 

( 1.  0.  2.  0.  0.) 

( 0,  3,  0.  0.  0.) 

M - ( 0.  4.  5.  6.  0.) 

( 0.  0.  0.  7.  0.) 


I 


c 

would  be  scored  as 

c 

1 

1 2 

3 

4 

5 6 

c 

— 

— , 

c 

lA  1 

1 3 

4 

7 

d 10 

c 

JA  1 

1 3 

2 

2 

3 4 

c 

A 1 

1.  2. 

3. 

4. 

S.  6 

The  strict  triangular  portions  of  the  matrices  L and  U are 
scored  in  Che  same  fashion  using  the  arrays  IL,  JL,  L and 
lU,  JU,  U respectively.  The  diagonal  entries  of  L and  U are 
assumed  to  be  equal  to  one  and  are  not  stored.  The  array  D 
contains  Che  reciprocals  of  Che  diagonal  entries  of  Che  matrix  D. 

IV.  ADDITIONAL  STORAGE  SAVINGS 

In  NSF,  R and  IC  can  be  the  same  array  in  the  calling 
sequence  if  no  reordering  of  the  coefficient  matrix  has  been  done. 

In  NNF,  R,  C and  1C  can  all  be  the  same  array  if  no  reordering 
has  been  done.  If  only  the  rows  have  been  reordered,  then  C and  IC 
can  be  the  same  array.  If  Che  row  and  column  orderings  are  the 
same,  then  R and  C can  be  Che  same  array.  Z and  ROW  can  be  the 
same  array. 

In  NNS,  R and  C can  be  the  same  array  if  no  reordering  has 
been  done  or  if  the  row  and  column  orderings  are  the  same.  Z and  B 
can  be  the  same  array;  however,  then  B will  be  destroyed. 

V.  PARAMETERS 

Following  is  a list  of  parameters  to  the  programs.  Names  are 
uniform  among  Che  various  subroutines.  Class  abbreviations  are  -- 
n - INTEGER  variable 
f - REAL  variable 

v - supplies  a VALUE  to  a subroutine 
r - returns  a RESULT  from  a subroutine 
i - used  INTERNALly  by  a subroutine 
a - ARRAY 

Class  I Parameter 


fva 


fva 


fvra 


FLAG 


- nonzero  entries  of  the  coefficient  matrix  M,  stored 

by  rows. 

Size  ■ number  of  nonzero  entries  in  M. 

- right-hand  side  b. 

Size  - N. 

- ordering  of  Che  columns  of  M. 

Size  • N. 

- reciprocals  of  the  diagonal  entries  of  the  matrix  D. 

Size  - N. 

- error  flag;  values  and  their  meanings  are  — 

0 No  Errors  Detected 

N-MC  Null  Row  in  A ~ Row  - K 

Duplicate  Entry  in  A — Row  - K 


2N+R 

3N-HC 

4N+1 

5N+K 

6N+K 

7N+1 

8N-HC 


lA 


IC 


Insufficient  Storage  for  JL  — Row  - K 
Insufficient  Storage  for  L 
Null  Pivot  — Row  « K 
Insufficient  Storage  for  JU  — Row  ■ K 
Insufficient  Storage  for  U 
Zero  Pivot  — Row  • K 

- pointers  to  delimit  the  rows  in  A. 

Size  - N+l. 

- Inverse  of  the  ordering  of  Che  columns  of  M;  i.e. 

1C(C(1)  - 1 for  I-l, ...N. 

Size  - N. 

- pointers  to  delimit  the  rows  in  L. 

Size  - N+l. 


IL 


C nvra 
C 

C nva 
C 

C nvra 
C 

C nv 

C 

C 

C nvra 


C nva 
C 

C fvra 


pointers  to  dellnit  the  rows  in  U, 

Sire  - N+l. 

column  numbers  corresponding  to  the  elements  oi  A, 

Size  • size  of  A. 

column  numbers  corresponding  to  the  elements  of  L. 

Size  - JLMAX. 

declared  dimension  of  JL;  JLMAX  must  be  larger  than 
the  number  of  nonzero  entries  in  the  strict  lower 
triangle  of  H plus  flllin  (ILlN+l)-!  after  NSF). 

column  numbers  corresponding  to  the  elements  of  U. 

Size  - JUHAX. 

declared  dimension  of  JU ; JUHAX  ssist  be  larger  than 
the  number  of  nonzero  entries  in  the  strict  upper 
triangle  of  H plus  flllin  (lU(Nel)-l  after  NSF). 

nonzero  entries  In  the  strict  lower  triangular  portion 
of  the  matrix  L,  stored  by  rows. 

Size  - LKAX 

declared  dimension  of  L;  LHAX  must  be  larger  than 
the  number  of  nonzero  entries  In  the  strict  lower 
triangle  of  H plus  flllin  (lL(N‘fl)-l  after  NSF). 

nisaber  of  varlables/equatlons. 

ordering  of  the  rows  of  M. 

Size  - N. 

nonzero  entries  In  the  strict  upper  triangular  portion 
of  the  matrix  U,  stored  by  rows. 

Size  - UMAX. 

declared  dimension  of  U;  UMAX  must  be  larger  t)ian 
the  number  of  nonzero  entries  In  the  strict  upper 
triangle  of  M plus  flllin  (lU(N-Fl)-l  after  NSF). 

solution  X. 

Size  - N. 


C***  Subroutine  NSF 

C***  Symbolic  LDU-factorlzatlon  of  a nonsymmetrlc  sparse  matrix 
C (uncompressed  pointer  storage) 

C 

SUBROUTINE  NSF 

* (N,  R,IC,  lA.JA,  IL,JL, JLMAX,  lU.JU.JUMAX,  Q,  IM,  FLAG) 


Input  variables: 
Output  variables: 


N.  R,  IC,  LA.JA,  JLMAX,  JUMAX. 
IL,JL,  1U,JU,  FLAG, 


Parameters  used  Internally: 

I Q - suppose  M'  Is  the  result  of  reordering  H;  If 
I processing  of  the  Kth  row  of  M'  (hence  the  Kth  rows 

I of  L and  U)  Is  being  done,  then  Q(J)  Is  Initially 

I nonzero  If  N'(K,J)  Is  nonzero;  since  values  need 

I not  be  stored,  each  entry  points  to  the  next 

I nonzero;  for  example.  If  N-9  and  the  3th  row  of 

M'  Is 


OxxOxOOxO, 


1 


nla 


IM 


then  Q will  Initially  be 

aJ^a8aal0a2  (a-  arbitrary); 

Q(N-ft ) points  to  the  first  nonzero  in  the  row  and 
the  last  nonzero  points  to  N+l ; as  the  algorithm 
proceeds,  other  elements  of  Q are  inserted  in  the 
list  because  of  fillln. 

Size  - N+1. 

in  the  factorization,  IM(1)  is  the  last 
the  Ith  row  of  U which  needs  to  be 
in  computing  fill  in. 


at  each  step 
element  in 
considered 
Size  - N. 


Internal  variables — 

JLPTR  - points  to  the  last  position  used  in  JL. 
JUPTR  - paints  to  the  last  position  used  in  JU. 


INTEGER  R(l),  IC(1),  IA(1).  JA(1),  1L(1),  JL(1), 

* 1U(1),  JU(1),  Qd),  IM(1),  FLAG,  (Jt,  VJ 

C 

******  In  it idl iz6  po in t6r6  **************************************** 

JLPTR  - 0 
IL(1)  - 1 
JUPTR  - 0 
IU(1)  - 1 
C 

******  For  each  row  of  L and  U ************************************ 

DO  10  K-1,N 

Q ******  Set  Q to  the  reordered  row  of  A **************************** 

Q(N+1)  - N+1 
JMIN  - IA(R(K)) 

JMAX  - 1A(R(K)+1 ) - 1 
IF  (JMIN. GT. JMAX)  GO  TO  101 
DO  2 J-JMIN.JMAX 
VJ  - 1C(JA(J)) 

QM  - N + 1 

1 M - QM 

gw  - Q(M) 

IF  (QM.LT.VJ)  GO  TO  1 
IF  (QM.EQ.VJ)  GO  TO  102 
Q(M)  - VJ 
Q(VJ)  - QM 

2 CONTINUE 


C 

(;  For  each  entry  in  the  lower  triangle  ***************+***++*+ 

1 - N+1 

3 I - Q(I) 

IF  (I.GE.K)  GO  TO  7 

C ******  l(K,I)  will  be  nonzero,  so  add  it  to  JL  ******************** 
JLPTR  - JLPTR+1 

IF  (JLPTR, GT.JLMAX)  GO  TC  03 
JL(JLPTR)  - I 

- I 


! 


t 


C ******  Inspect  Ith  row  for  fillln,  adjust  IM  if  possible  ********** 
JMIN  - 1U(I) 

JMAX  - IM(I) 

IF  (JMIN. CT. JMAX)  GO  TO  6 
DO  5 J*JMIN,JMAX 
VJ  - JU(J) 

IF  (VJ.EQ.K)  1M(I)  - J 

4 M - QM 

QM  - Q(M) 

IF  (QM.LT.VJ)  GO  TO  4 
IF  (QM.EQ.VJ)  GO  TO  5 
Q(M)  - VJ 
q(VJ)  - QM 
QM  - VJ 

5 CONTINUE 

6 GO  TO  3 
C 

******  Check  for  null  pivot  *************************************** 

7 IF  (I.NE.K)  GO  TO  105 

C ******  Remaining  elements  of  Q define  structure  of  U(K,  ) ********* 

8 I - Q(I) 

IF  (l.GT.N)  GO  TO  9 

JUPTR  - JUPTR+1 

IF  (JUPTR. GT.JUMAX)  GO  TO  106 
JU  (JUPTR)  - I 
GO  TO  8 

Q ******  Get  reedy  for  next  row  ************************************* 

9 IK(K)  - JCfPTR 
IL(K+1)  - JLPTR+l 

10  IU(K+1)  - JUPTR+1 

C 

FLAG  - 0 
RETURN 
C 

C **  ERROR:  Null  Row  in  A 

101  FLAG  - N + R(K) 

RETURN 

C **  ERROR:  Duplicate  Entry  in  A 

102  FLAG  - 2*N  + R(K) 

RETURN 

C **  ERROR:  Insufficient  Storage  for  JL 

103  FLAG  - 3*N  + K 
RETURN 

C **  ERROR:  Null  Pivot 

105  FLAG  - 5*N  + K 
RETURN 

C **  ERROR:  Insufficient  Storage  for  JU 

106  FLAG  - 6*N  + K 
RETURN 


c 

c 

c 

C***  Subroutine  NNF 
C***  Numeric  LDU-factorizatlon  of  sparse  nonsymmecrlc  matrix  and 
C solution  of  system  of  linear  equations  (uncompressed  pointer 

C storage) 

C 

SUBROUTINE  NNF 

* (N,  R.C.IC,  U.JA.A,  Z,  B,  IL, JL, L, LMAX,  D,  lU, JU, U, UMAX, 

* ROW,  IMP,  FLAG) 

C 

C Input  variables:  N,  R,C,IC,  LA,JA,A,  B,  IL,JL,LMAX,  1U,JU,UMAX 

C Output  variables:  Z,  L,D,U,  FLAG 

C 

C Parameters  used  internally: 

C fla  I ROW  - holds  Intermediate  values  In  calculation  of  L,  D,  U. 

C I Size  » N. 

C fla  I TMP  - holds  new  right-hand  side  b'  for  solution  of  the 
C I equation  Ux  “ b' . 

C I Size  ” N. 

C 

INTEGER  R(l),  C(l),  IC(1),  IA(1),  JA(1), 

* IL(l),  JL(1),  LMAX,  IU(1),  JU(1),  UMAX,  FUG 

REAL  AO),  ZO),  BO),  LO),  DO),  UO),  ROWO),  TMPO),  LI 
C 

C ******  Check  storage  ********************************************** 

IF  (IL(N+1)-1  .GT.  LMAX)  GO  TO  104 
IF  (1U(N+1)-1  .GT.  UMAX)  GO  TO  107 
C 

C ******  For  each  row  *********************************************** 

DO  10  K-1,N 

C ******  Set  the  Initial  structure  of  ROW  *************************** 

JMIN  - IL(K) 

UMAX  - IL(K+1)  - 1 
IF  (JMIN. GT. UMAX)  GO  TO  2 

C ******  If  l(K,M)  .NE.  0,  ROW(M)-0  ********************************* 

DO  1 J-JMIN,JMAX 

1 R0W(JL(J))  - 0 

2 ROW(K)  - 0 
JMIN  - IU(K) 

JMAX  - IU(K+1)  - 1 
IP  (JMIN.  GT.  JMAX)  GO  TO  4 

C ******  F£  U(K,M)  .NE.  0,  R0W(M)-0  ********************************* 

DO  3 J-JMIN,JMAX 

3 R0W(JU(J))  - 0 

4 JMIN  - IA(R(K)) 

JMAX  - IA(R(K)+l)  - 1 

C ******  Set  ROW  to  Kth  row  of  reordered  A ************************** 

DO  5 J«JMIN,JMAX 

5 R0W(1C(JA(J)))  - A(J) 

C ******  Initialize  SUM  ********************************************* 

SUM  - B(R(K)) 


non  nn  no  n nn 


******  Assign  the  Kth  row  of  L and  adjust  ROW,  SUM  **************** 
IMIN  - IL(K) 

IMAX  - IL(K+1)  - 1 
IF  (IMIN. GT. IMAX)  GO  TO  8 
DO  7 I-IMIN,IMAX 
LI  - - ROW(JL(I)) 

******  If  L Is  not  required,  then  comment  out  the  following  line  ** 
L(l)  - - LI 

SUM  - SUM  + LI  * TMP(JL(1)) 

JMIN  - IU(JL(I)) 

JMAX  - IU(JL(I)+1)  - 1 
IF  ( JMIN. GT. JMAX)  GO  TO  7 
DO  6 J"JMIN,JMAX 

6 ROW(JU(J))  - ROW(JU(J))  + LI  * 0(J) 

7 CONTINUE 

******  Assign  diagonal  D and  Kth  row  of  U,  set  TMP(K)  ************* 

8 IF  (ROW(K).EQ.O)  GO  TO  108 
DK  - 1 / ROW(K) 

D(K)  - DK 
TMP(K)  - SUM  * DK 
JMIN  - 1U(K) 

JMAX  - 1U(K+1)  - I 
IF  (JMIN. GT. JMAX)  CO  TO  10 
DO  9 J-JMIN,JMAX 

9 U(J)  - R0W(JU(J))  * DK 

10  CONTINUE 

******  Solve  Ux  - TMP  by  back  substitution  ********************** 
K - N 

DO  13  1-1, N 
SUM  - TMP(K) 

JMIN  - 1U(K) 

JMAX  - IU(K+1)  - 1 
IF  (JMIN. GT. JMAX)  GO  TO  12 
DO  11  J-JMIN,JMAX 

11  SUM  - SUM  - U(J)  * TMP(JU(J)) 

12  TMP(K)  - SUM 
Z(C(K))  - SUM 

13  K - K-1 

FLAG  - 0 
RETURN 

**  ERROR:  Insufficient  Storage  for  L 
104  FLAG  - 4*N  + 1 
RETURN 

C **  ERROR:  Insufficient  Storage  for  U 

107  FLAG  - 7*N  + 1 
RETURN 

C **  ERROR:  Zero  Pivot 

108  FLAG  - 8*N  + K 
RETURN 

END 


onooonnno  nnooonn 


***  Subroutine  NNS 

k**  Numeric  solution  of  a sparse  nonsymmetrlc  system  of  linear 

equations  given  LDU-factorizatlon  (uncompressed  pointer  storage) 

SUBROUTINE  NNS 

* (N.  R,C.  IL.JL.L,  D,  lU.JU.U,  Z,  B,  TMP) 

Input  variables:  N,  R,C,  IL.JL.L,  D,  lU.JU.U,  B 

Output  variables:  Z 

Parameters  used  Internally: 

f la  I TMP  - holds  new  right-hand  side  b'  for  solution  of  Che 

I equation  Ux  “ b' . 

I Size  - N. 

INTEGER  R(l),  C(l),  1L(1),  JL(1),  IU(1),  JU(1) 

REAL  Ld).  D(l).  U(l),  Z(l),  B(l),  TMP(l) 

******  Solve  LDy  * b by  forward  substitution  ********************* 

DO  2 K-1,N 
SUM  - B(R(K)) 

JMIN  - IL(K) 

JMAX  - 1L(K+1)  - 1 
IF  (JMIN. GT. JMAX)  GO  TO  2 
DO  I J-JMIN,JMAX 

1 SUM  - SUM  - L(J)  * TMP(JL(J)) 

2 TMP(K)  "SUM  * D(K) 

C 

C ******  Solve  Ux  » y by  back  substitution  ************************ 

K - N 

DO  5 1-1,  N 
SUM  - TMP(K) 

JMIN  - IJ(K) 

JMAX  - IU(K+1)  - I 
IF  (JMIN. GT.  JMAX)  GO  TO  4 
DO  3 J-JMIN,JMAX 

3 SUM  - SUM  - U(J)  * TMP(JU(J)) 

4 TMP(K)  - SUM 
Z(C(K))  - SUM 

5 K - K-l 

RETURN 


Appendix  2 
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C 
C 

C Subroutines  for  Solving  Sparse  Nonsynmetrlc  Systems 

C of  Linear  Equations  (Track  Nonzeroes  Dynamically) 

C 

C 

C***  Subroutine  TDRV 

C***  Driver  for  subroutine  for  solving  sparse  nonsymmetrlc  systems  of 
C linear  equations  (track  nonzeroes  dynamically) 

C 

SUBROUTINE  TDRV 

* (N.  R,IC,  lA.JA.A,  B,  Z.  NSP, ISP, RSP, ESP,  FLAG) 

C 

C PARAMETERS 

C Class  abbreviations  are  — 

C n - INTEGER  variable 

C f - r.EAL  variable 

C V - supplies  a VALUE  to  the  driver 

C r - returns  a RESULT  from  the  driver 

C 1 - used  INTERNALly  by  the  driver 

C a - ARRAY 

C 

C Class  I Parameter 

C ( 

C I 

C The  nonzero  entries  of  the  coefficient  matrix  M are  stored 

C row-by-row  In  the  array  A.  To  Identify  the  Individual  nonzero 
C entries  In  each  row,  we  need  to  know  in  which  column  each  entry 

C lies.  The  column  Indices  which  correspond  to  the  nonzero  entries 

C of  M are  stored  In  the  array  JA;  l.e.,  if  A(K)  - H(1,J),  then 
C JA(K)  ■ J.  In  addition,  we  need  to  know  where  each  row  starts  and 

C how  long  It  Is.  The  Index  positions  in  JA  and  A where  Che  rows  of 

C N begin  are  stored  In  the  array  LA;  l.e.,  if  M(1,J)  Is  the  first 

C nonzero  entry  (stored)  in  the  I-th  row  and  A(K)  - M(I,J),  then 

C IA(I)  " K.  Moreover,  the  Index  In  JA  and  A of  the  first  location 

C following  the  last  element  In  the  last  row  Is  stored  In  IA(N-I-1 ) . 

C Thus,  the  number  of  entries  In  the  1-th  row  is  given  by 

C IA(1+1)  - IA(1),  the  nonzero  entries  of  the  1-th  row  are  stored 

C consecutively  In 

C A(1A(I)),  A(IA(l)+l),  ....  A(IA(I+1)-1), 

C and  Che  corresponding  column  Indices  are  stored  consecutively  In 
C JA(IA(I)).  JA(1A(I)+1),  ....  JA(1A(I+1)-1). 

C For  example,  the  5 by  3 matrix 


C 

( 

1.  0.  2.  0.  0.) 

c 

( 

0.  3.  0.  0.  0.) 

c 

M - ( 

0.  4.  5.  6.  0.) 

c 

( 

0.  0.  0.  7.  0.) 

c 

( 

0.  0.  0.  8.  9.) 

C would  be  stored  as 

C 1123456789 

C 1 

C lA  I 1 3 4 7 8 10 

C JA  1132234445 

C A I 1.  2,  3.  4.  5,  6.  7.  8.  9. 

C 

C nv  I N - number  of  variables/ equations. 

C fva  I A - nonzero  entries  of  the  coefficient  matrix  M,  stored 

C I by  rows . 

C I Size  • niaiber  of  nonzero  entries  In  M. 

C nva  I lA  - pointers  to  delimit  the  rows  In  A. 


r 


c 

1 

Size  - N+1. 

c 

nva 

1 JA 

- column  numbers  corresponding  to 

the  elements  of 

c 

1 

Size  - size  of  A. 

c 

f va 

1 B 

- right-hand  side  b;  B and  Z can 

Che  same  array. 

c 

1 

Size  - N. 

c 

fra 

1 z 

- solution  x;  B and  Z can  be  Che 

same  array. 

c 

1 

Size  • N. 

The  rows  and  columns  of  the  original  matrix  M can  be 
reordered  (e*g.,  to  reduce  illlin  or  ensure  numerical  stability) 
before  calling  Che  driver.  If  no  reordering  is  done,  then  set 
R(I)  • C(I)  • IC(I)  - I for  1*1, The  solution  Z is  returned 
in  Che  original  order. 


- ordering  of  the  rows  of  M. 


Inverse  of  the  ordering  of  the  columns  of  H; 


IC(C(I))  - I for  I-l N, 

ordering  of  the  columns  of  M. 
Size  - N. 


where  C is  the 


Various  errors  are  detected  by  the  driver  and  the  individual 
subroutines. 


FLAG  - error  flag;  values  and  their  meanings  are  -- 
0 No  Errors  Detected 
N-rtC  Null  Row  in  A ~ Row  - K 
2N+IC  Duplicate  Entry  in  A — Row  - K 
5N+K  Null  Fivot  --  Row  - K 

8N+K  Zero  Pivot  --  Row  ■ K 

lON+1  Insufficient  Storage  in  TDRV 

IZN-fK  Insufficient  Storage  in  TRK 


Working  storage  is  needed  for  the  factored  form  of  the  matrix 
M plus  various  temporary  vectors.  The  arrays  ISP  and  RSP  should  be 
the  same;  integer  storage  is  allocated  from  the  beginning  of  ISP 
and  real  storage  from  the  end  of  RSP. 


C 

C 

C nvira 
C 
C 
C 

C fvlra 


declared  dimension  of  ISP  and  RSP;  NSP  generally  must 
be  larger  than  6N-t-2  + 2*K  (where  K - (nisaber  of 
nonzero  entries  In  the  upper  triangle  of  M) ) . 

Integer  working  storage  divided  up  into  various  arrays 
needed  by  the  subroutines;  ISP  and  RSP  should  be 
the  same  array. 

Size  - NSP. 

real  working  storage  divided  up  into  various  arrays 
needed  by  the  subroutines;  ISP  and  RSP  should  be 
the  same  array. 

Size  - NSP. 

if  NSP  is  sufficiently  large  to  allocate  space,  then 
ESP  is  set  to  the  amount  of  excess  storage  provided. 


INTEGER  R(l),  IC(1),  1A(1),  JA( 

U,  ROW,  TUP,  Q 

REAL  A(l),  B(l),  Z(l),  RSP(l) 


1A(1),  JA(l),  ISP(l),  ESP, 


w 


***************** 


c 

C ******  Initialize  and  divide  up  temporary  storage 


IJU 

1 

lU 

IJU 

+ N 

Q 

lU 

+ N+1 

IM 

Q 

+ N+1 

JU 

IM 

+ N 

U 

JU 

ROW 

NSP 

- N 

TMP 

ROW 

- N 

MAX 

TMP 

- JU 

IF 

(MAX.LT.O)  GO  TO  110 

C ******  Call  zero-tracking  subroutine  ****************************** 

FLAG  - 0 
CALL  TRK 

* (N,  R,  IC,  lA,  JA.  A,  Z,  B, 

* ISP(IJU),  ISP(JU).  ISP(IU),  RSP(U),  MAX, 

* ISP(Q),  ISP(IM),  RSP(ROW),  RSP(TMP),  FLAG,  ESP) 

IF  (FLAG.NE.O)  GO  TO  100 

RETURN 

C 

C **  ERROR:  Error  Detected  In  TRK 
100  RETURN 

C **  ERROR:  Xnsulflclent  Storage 
110  FLAG  • 10»N  + 1 
RETURN 
END 


YALE  SPARSE  MATRIX  PACKAGE  - ZERO-TRACKING  CODE 
SOLVING  THE  SYSTEM  OF  EQUATIONS  Mx  - b 
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0 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


I.  SUBROUTINE  NAMES 

TRK  performs  an  LDU-decomposition  of  the  matrix  M,  without 
storing  L or  D,  and  solves  the  linear  system  of  equations. 


II.  CALLING  SEQUENCES 

The  coefficient  matrix  can  be  processed  by  an  ordering  routine 
(e.g.,  to  reduce  flllln  or  ensure  numerical  stability)  before  using 
Che  remaining  subroutines.  If  no  reordering  is  done,  then  set 
R(I)  ■ C(I)  ■ IC(I)  • I for  I-1,...,N.  The  calling  sequence  is  — 
( (matrix  ordering)) 

TRK  (solution  of  linear  system  of  equations) 

(If  several  systems  with  Che  same  coefficient  matrix  but  different 
right-hand  sides  or  several  systems  whose  coefficient  matrices  have 
Che  same  nonzero  structure  are  to  be  solved,  and  sufficient  space 
is  available,  other  subroutines  should  be  used.) 


III.  STORAGE  OF  SPARSE  MATRICES 

The  nonzero  entries  of  the  coefficient  matrix  M are  stored 
row-by-row  in  the  array  A.  To  identify  the  individual  nonzero 
entries  in  each  row,  we  need  to  know  in  which  column  each  entry 
lies.  The  column  Indices  which  correspond  to  the  nonzero  entries 
of  M are  stored  in  the  array  JA;  i.e.,  if  A(K)  ■ M(1,J),  then 
JA(K)  • J.  In  addition,  we  need  to  know  where  each  row  starts  and 
how  long  it  is.  The  index  positions  in  JA  and  A where  the  rows  of 
M begin  are  stored  in  the  array  lA;  i.e.,  if  M(I,J)  is  the  first 
nonzero  entry  (stored)  in  the  1-th  row  and  A(K)  « M(I,J),  then 
IA(I)  ■■  K.  Moreover,  the  index  in  JA  and  A of  the  first  location 
following  Che  last  element  in  the  last  row  is  stored  in  IA(N+1). 
Thus,  the  number  of  entries  in  the  I-th  row  is  given  by 
IA(I+1)  - IA(I),  the  nonzero  entries  of  the  I-th  row  are  stored 
consecutively  in 

A(IA(I)),  A(IA(1)+1),  ....  A(IA(I+1)-1), 

and  the  corresponding  column  indices  are  stored  consecutively  in 

JA(IA(1)),  JA(IA(I)+1),  ...,  JA(1A(I+1)-1). 

For  example,  the  3 by  3 matrix 

( 1.  0.  2.  0.  0.) 

( 0.  3.  0.  0.  0.) 

M - ( 0.  4.  3.  6.  0.) 

( 0.  0.  0.  7.  0.) 

( 0.  0.  0.  8.  9.) 
would  be  stored  as 

112  3 4 3 6 7 8 9 


lA  I I 3 4 7 8 10 

JA  1132234445 
A I 1.  2.  3.  4.  5.  6.  7.  8.  9. 

The  strict  upper  triangular  portion  of  the  matrix  U is  stored 
in  a similar  fashion  using  the  arrays  lU,  JU,  U,  except  that  an 
additional  array  IJU  is  used  to  compress  storage  of  JU  by  allowing 
some  of  the  column  indices  to  be  used  for  more  than  one  row. 

IJU(K)  points  to  the  starting  location  in  JU  of  entries  for  the  Kth 
row.  Compression  in  JU  occurs  in  two  ways.  First,  if  a row  I was 
merged  into  the  current  row  K,  and  the  number  of  elements  merged  in 
from  (the  tail  portion  of)  row  I is  the  same  as  the  final  length  of 
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i 


1 


t 


i 

I 

i 

J 

r 

j 

I 

1 

I 
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c 

row  K,  Chen  the  Kth  row  and  Che  Call  of  cow  1 are  Identical  and 

c 

IJU(K)  may  point  to  the  start  of  the  tall.  Second,  If  some  tall 

c 

portion  of  the  (K-1 )8t  row  Is  Identical  to  the  head  of  the  Kth  cow. 

c 

then  IJU(K)  may  point  to  the  start  of  that  tail  portion.  For 

c 

example,  Che  nonzero  structure  of  Che  strict  upper  triangular  part 

c 

of  the  matrix 

c 

d 0 X X X 

c 

0 d 0 X X 

c 

0 0 d X 0 

c 

0 0 0 d X 

c 

0 0 0 0 d 

c 

would  be  represented  as 

c 

11  2 3 4 5 6 

c 

lU  1 1 4 6 7 8 8 

c 

JU  13  4 5 4 

c 

IJU  1 1 2 4 3 

c 

IV.  ADDITIONAL  STORAGE  SAVINGS 

c 

JU  and  U should  be  Che  same  array.  TRK  fills  JU  from  the 

c 

beginning  of  Che  array  and  U from  the  end  of  the  array. 

c 

R and  1C  can  be  Che  same  array  in  the  calling  sequence  if  no 

c 

reordering  of  Che  coefficient  matrix  has  been  done.  Z and  ROW  can 

c 

be  Che  same  array. 

c 

V.  PARAMETERS 

c 

Following  Is  a list  of  parameters  Co  TRK.  Class  abbreviations 

c 

are  — 

c 

n - INTEGER  variable 

c 

f - REAL  variable 

c 

V - supplies  a VALUE  to  a subroutine 

c 

r - returns  a RESULT  from  a subroutine 

c 

1 - used  INTERNALly  by  a subroutine 

c 

a - ARRAY 

c 

Class  1 Parameter 

c 

» 

c 

Eva  1 A - nonzero  entries  of  the  coefficient  matrix  M,  stored 

c 

1 by  rows. 

c 

1 Size  • number  of  nonzero  entries  in  M. 

c 

fva  1 B - right-hand  side  b. 

c 

1 Size  » N. 

c 

nr  1 ESP  - if  enough  storage  was  provided  for  JU  and  U,  Chen  ESP 

c 

1 Is  set  to  amount  of  excess  storage  provided. 

c 

nr  1 FLAG  - error  flag;  values  and  their  meanings  are  — 

c 

1 0 No  Errors  Detected 

c 

1 N-4C  Null  Row  in  A — Row  - K 

c 

1 2N-tK  Duplicate  Entry  in  A — Row  - K 

c 

1 5N+K  Hull  Pivot  — Row  » K 

c 

1 8N-HC  Zero  Pivot  — Row  - K 

c 

1 12N-fK  Insufficient  Storage  for  JU/U  — Row  • K 

c 

nva  1 lA  - pointers  to  delimit  the  rows  In  A. 

c 

1 Size  " N+1. 

c 

nva  1 IC  - Inverse  of  the  ordering  of  the  columns  of  M;  l.e.. 

c 

1 IC(C(1))  - 1 for  I-l,...N,  where  C is  the 

. j 

_ 

r ' 


ordering  of  the  columns  of  M. 

Size  - N. 

- pointers  to  the  first  element  in  each  row  In  JU,  used 

to  compress  storage  In  JU. 

Size  - N. 

- pointers  to  delimit  the  rows  in  U. 

Size  - N+1. 

- column  numbers  corresponding  to  the  elements  of  A. 

Size  - size  of  A. 

- column  numbers  corresponding  to  the  elements  of  U; 

JU  and  U should  be  the  same  array. 

Size  - MAX. 

- declared  dimension  of  JU  and  U;  MAX  must  be  larger 

chan  Che  size  of  U (the  number  of  nonzero  entries 
in  the  strict  upper  triangle  of  M plus  fillln)  plus 
Che  size  of  JU  (the  size  of  U minus  compression). 

- number  of  variables/equations. 

- ordering  of  the  rows  of  M. 

Size  - N. 

- nonzero  entries  In  the  strict  upper  triangular  portion 

of  U,  stored  by  rows;  JU  and  U should  be  the  same 
array. 

Size  - MAX. 

- solution  X. 

Size  » N. 


C***  Subroutine  TRK 

C***  Numerical  solution  of  sparse  nonsymmetric  system  of  linear 
C equations  (track  zeroes  dynamically) 

C 

SUBROUTINE  TRK 

* (N,  R,IC,  lA.JA.A,  Z,  B.  IJU, JU, lU, U.MAX, 

* Q,  IM,  ROW,  TMP,  FLAG,  ESP) 

C 

C Input  variables;  N,  R, IC,  IA,JA,A,  B,  MAX 

C Output  variables:  Z,  FLAG 


Parameters  used  Internally: 

I Q - suppose  M'  is  the  result  of  reordering  M;  if 
I processing  of  the  Kth  row  of  M'  (hence  the  Kth  rows 

I of  L and  U)  is  being  done,  then  Q(J)  is  initially 

I nonzero  if  M'(K,J)  is  nonzero;  since  values  need 

I not  be  stored,  each  entry  points  to  the  next 

I nonzero;  for  example,  if  N-9  and  the  5th  row  of 

I M'  is 

I OxxOxOOxO, 

I then  Q will  initially  be 

I a35a8aal0a2  (a-  arbitrary); 

I Q(N-Pl)  points  to  the  first  nonzero  in  the  row  and 

I the  last  nonzero  points  to  N-fl ; as  the  algorithm 

I proceeds,  other  elements  of  Q are  Inserted  in  the 

I list  because  of  fillln. 

I Size  - N+1. 

I IM  - at  each  step  in  the  factorization,  IM(I)  is  Che  last 
I element  of  the  Ith  row  of  U which  needs  to  be 

I considered  in  computing  fillln. 

Size  - N. 


d 


t 


1 


f la 


f la 


ROW  - holds  Intermediate  values  In  calculation  of  U. 
Size  “ N. 

IMP  - holds  new  right-hand  side  b'  for  solution  of  the 
equation  Ux  - b' . 

Size  • N. 


INTEGER  R(l),  IC(1),  1A(1).  JAd), 

* IJU(l),  JU(l).  IU(1),  Q(l),  IM(1),  FLAG,  ESP,  VJ, 

REAL  A(l),  Z(l),  B{1),  U(l),  ROU(l),  TMP(l) 

C 

C ******  Initialize  ************************************************* 


JUMIN  - 1 
JUMAX  - 0 
IU(1)  - MAX 
C 

******  each  row  a********************************************** 


DO  20  K-1,N 

C ******  Initialize  Q and  ROW  to  the  Xth  row  of  reordered  A ********* 
LUX  - 0 
Q(N+1)  - N+1 
JMIN  - 1A(R(X)) 

JMAX  - 1A(R(K)+1)  - 1 
IF  (JMIN. GT. JMAX)  GO  TO  101 
DO  2 J-JMIN,JMAX 
VJ  - IC(JA(J)) 
qM  - N+1 

1 M - QM 
OH  - q(M) 

IF  (QM.LT.VJ)  GO  TO  1 
IF  (QM.EQ.VJ)  GO  TO  102 
LUK  • LUK+1 
Q(M)  = VJ 
Q(VJ)  - QM 
ROWCVJ)  - A(J) 

2 CONTINUE 


C 

Q ******  Link  through  Q ********************************************* 
LHAX  - 0 
IJU(K)  - JUMAX 
I - N+1 

3 I - Q(I) 

LUK  - LUK-1 

IF  (I.GE.K)  GO  TO  8 
OM  - I 

JMIN  - IJU(I) 

JMAX  - IM(I) 

LUI  - 0 

IF  (JMIN. GT. JMAX)  GO  TO  7 

C ******  and  find  nonzero  structure  of  Kth  row  of  L and  U *********** 
DO  3 J>JM1N,JMAX 
VJ  - JU(J) 

IF  (VJ.GT.K)  LUI  - LUI+1 

4 M - QM 

qn  - Q(M) 

IF  (QM.LT.VJ)  GO  TO  4 
IF  (QM.EQ.VJ)  GO  TO  5 
LUK  - LUK+1 
Q(M)  - VJ 
Q(VJ)  - QM 
ROW(VJ)  - 0 
QM  - VJ 

5 CONTINUE 


1 

i 

■J 

i 


r ' 


Q ******  Adjust  IJU  and  IM  ****************************************** 

JTMP  - JMAX  - LUI 
IF  (LUl.LE.LMAX)  GO  TO  6 
LMAX  - LUI 
IJU(IC)  - JTMP+1 

6 IF  (JTMP.LT.JMIN)  GO  TO  7 

IF  (JU(JTMP).EQ.K)  IM(I)  - JTMP 

7 GO  TO  3 


C 

C ******  See  if  JU  storage  can  be  compressed  ************************ 

8 IF  (I.NE.K)  GO  TO  105 

IF  (LUK.EQ.LMAX)  GO  TO  14 
IF  (JUMIN.GT.  JUMAX)  GO  TO  12 
1 - Q(K) 

DO  9 JMIN-JUMIN,  JUMAX 

IF  (JU(JMIN)-I)  9,  10,  12 

9 CONTINUE 
GO  TO  12 

10  IJU(K)  - JMIN 

DO  11  J-JMIN, JUMAX 

IF  (JU(J).NE.I)  GO  TO  12 
1 - Q(I) 

IF  (I.GT.N)  GO  TO  14 

11  CONTINUE 
JUMAX  - JMIN  - 1 

Q ******  Store  pointers  in  JU  *************************************** 

12  JUMIN  ■=  JUMAX  + 1 

JUMAX  - JUMAX  + LUK 

IF  ( JUMAX. GT.IU(K))  GO  TO  112 
I - K 

DO  13  J-JUMIN,  JUMAX 

I = Q(l) 


13 

JU(J)  - I 

IJU(K)  - JUMIN 

14 

IU(K+1)  - IU(K)  - LUK 

IF  (JUMAX.GT.IU(K+l))  GO  TO  112 

IM(K)  - IJU(K)  + LUK  - 1 

****** 

Calculate  numerical  values  for  Kth  row  ********************* 

SUM  - B(R(K)) 

I - N+1 

15 

I - Qd) 

IF  (I.GE.K)  GO  TO  18 

AKI  - - ROW(I) 

SUM  - SUM  + AKI  * TMP(I) 

JMIN  - IU(I+1)  + 1 

JMAX  - 1U(I) 

IF  (JMIN.  GT.  JMAX)  GO  TO  17 

MU  - IJU(I)  - JMIN 

16 

17 

DO  16  J-JM1N,JMAX 

R0«(JU(MU+J))  - ROW(JU(MU+J))  + AKI  * U(J) 

GO  TO  15 

16 

17 


C ******  score  values  In  TMP  and  U ********** 

18  IF  (ROW(K).EQ.O)  GO  TO  108 
DK  - 1 / ROW(IC) 

TMP(K)  - SUM  * OK 
JMIN  - 1U(K+1)  + 1 
JMAX  - IU(K) 

IF  (JMIN. GT. JMAX)  CO  TO  20 
MU  - IJU(K)  - JMIN 
DO  19  J-JMIN,JMAX 

19  U(J)  - ROW(JU(MU+J))  * DK 

20  CONTINUE 

ESP  - IU(N+1)  - JUMAX 
C 

C ******  Solve  Ux  “ TMP  by  back  substitution 
K - N 

DO  23  I-1,N 
SUM  - TMP(K) 

JMIN  - IU(K+1)  + 1 
JMAX  - IU(K) 

IF  ( JMIN. GT. JMAX)  GO  TO  22 
MU  - IJU(K)  - JMIN 
DO  21  J*JMIN,JMAX 

21  SUM  - SUM  - U(J)  * TMP(JU(MU+J)) 

22  TMP(K)  - SUM 

23  K - K-1 
DO  24  K-1,N 

24  Z(K)  - TMP(IC(K)) 

C 

FUG  - 0 
RETURN 
C 

C **  ERROR:  Null  Row  in  A 

101  FUG  - N + R(K) 

RETURN 

C **  ERROR:  Duplicate  Entry  in  A 

102  FUG  - 2*N  + R(K) 

RETURN 

C **  ERROR:  Null  Pivot 
105  FUG  - 5*N  + K 
RETURN 

C **  ERROR:  Zero  Pivot 
108  FUG  - 8*N  + K 
RETURN 


C **  ERROR:  Insufficient  Storage  for  JU  and  U 
112  FUG  - 12*N  + K 


I 


Appendix  i 


7/31/77 


Subroutines  for  Solving  Sparse  Nonsynnetrlc  Systems 
of  Linear  Equations  (Compressed  Pointer  Storage) 


C***  Subroutine  CDRV 

C**s  Driver  for  subroutines  for  solving  sparse  nonsymmetrlc  systems  of 
C linear  equations  (compressed  pointer  storage) 

C 

SUBROUTINE  CDRV 

* (N,  R,C,1C,  lA.JA.A,  B.  Z.  NSP, ISP, RSP, ESP,  PATH,  FLAG) 


C 

C PA 

C Cl 

C 

C 

C 

C 

C 

C 

C 

C Class 


PARAMETERS 

Class  abbreviations  are — 
n - INTEGER  variable 
f - REAL  variable 

V - supplies  a VALUE  to  Che  driver 
r - returns  a RESULT  from  Che  driver 
1 - used  INTERNALly  by  the  driver 
a - ARRAY 

iss  I Parameter 


The  nonzero  entries  of  the  coefficient  matrix  M are  scored 
row-by-row  In  the  array  A>  To  identify  the  individual  nonzero 
entries  in  each  row,  we  need  Co  know  in  which  column  each  entry 
lies.  The  column  Indices  which  correspond  to  Che  nonzero  entries 
of  M are  scored  in  the  array  JA;  i.e.,  if  A(K)  ~ M(1,J),  then 
JA(K)  - J.  In  addition,  we  need  to  know  where  each  row  starts  and 
how  long  it  is.  The  index  positions  in  JA  and  A where  the  rows  of 
M begin  are  scored  in  Che  array  lA;  i.e.,  if  M(1,J)  is  the  first 
nonzero  entry  (stored)  in  the  I-th  row  and  A(K)  ■■  M(I,J),  then 
IA(I)  • K.  Moreover,  the  index  in  JA  and  A of  the  first  location 
following  the  last  element  in  the  last  row  is  stored  in  lA(N-t-l). 
Thus,  Che  number  of  entries  in  the  1-th  row  is  given  by 
lA(I-t-l)  - 1A(I),  the  nonzero  entries  of  the  I-th  row  are  stored 
consecutively  In 

A(1A(I)),  A(IA(I)+1),  ....  A(IA(I+1)-1), 

and  the  corresponding  column  Indices  are  stored  consecutively  in 

JA(IA(I)),  JA(IA(I)+l),  JA(1A(I+1)-1). 

For  example,  the  5 by  3 matrix 

( 1.  0.  2.  0.  0.) 

( 0.  3.  0.  0.  0.) 

M - ( 0.  4.  5.  6,  0.) 

( 0.  0.  0.  7.  0.) 

( 0.  0.  0.  8.  9.) 
would  be  stored  as 

11  2 3 4 5 6 7 8 9 


1 3 4 7 8 10 

132234443 

1.  2.  3.  4.  5.  6,  7.  8.  9. 


number  of  variables/equations. 

nonzero  entries  of  the  coefficient  matrix  M,  stored 
by  rows. 

Size  • number  of  nonzero  entries  in  M. 
pointers  to  delimit  the  rows  in  A. 


JA 


C 

C nva 
C 

L i’ va 
C 

C f ra 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C nva 
C 

C nva 
C 

C nva 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

i: 

c 

c 

c 

c 

c 

c 

c 

C nv 

C 

C 

C 

C 

C 

C 

C 


B 

Z 


Size  - N+1. 

- colunn  nunbers  corresponding  to  Che  elements  of  A. 

Size  ••  size  of  A. 

- right-hand  side  b;  B and  Z can  Che  same  array. 

Size  • N. 

- solution  x;  B and  Z can  be  Che  same  array. 

Size  - N. 


The  rows  and  columns  of  Che  original  matrix  M can  be 
reordered  (e.g.,  Co  reduce  flllln  or  ensure  numerical  stability) 
before  calling  the  driver.  If  no  reordering  Is  done,  then  set 
R(l)  ~ C(I)  - IC(I)  - I for  1-1, ...,N.  The  solution  Z Is  returned 
In  the  original  order. 

If  Che  columns  have  been  reordered  (i.e.,  C(1).NE.I  for  some 

I),  Chen  the  driver  will  call  a subroutine  (NROC)  which  rearranges 
each  row  of  JA  and  A,  leaving  the  rows  in  Che  original  order,  but 
placing  Che  elements  of  each  row  In  Increasing  order  with  respect 
CO  Che  new  ordering.  If  PATH.NE. 1,  then  NROC  Is  assumed  to  have 
been  called  already. 

I R - ordering  of  the  rows  of  M. 

I Size  - N. 

1 C - ordering  of  the  columns  of  M. 

I Size  - N. 

I 1C  - Inverse  of  Che  ordering  of  the  columns  of  M;  l.e., 

I IC(C(I))  - I for  I-l N. 

j Size  - N. 

The  solution  of  the  system  of  linear  equations  is  divided  Into 
three  stages  — 

NSFC  — The  matrix  M Is  processed  symbolically  to  determine  where 
flllln  will  occur  during  Che  numeric  factorization. 

NNFC  — The  matrix  M Is  factored  numerically  into  Che  product  U>U 
of  a unit  lower  triangular  matrix  L,  a diagonal  matrix 
D,  and  a unit  upper  triangular  matrix  U,  and  the  system 
Mx  - b Is  solved. 

NNSC  — The  linear  system  Mx  - b is  solved  using  Che  LDU 
factorization  from  NNFC. 

For  several  systems  whose  coefficient  matrices  have  the  same 
nonzero  structure,  NSFC  need  be  done  only  once  (for  the  first 
system);  then  NNFC  Is  done  once  for  each  additional  system.  For 
several  systems  with  Che  same  coefficient  matrix,  NSFC  and  NNFC 
need  be  done  only  once  (for  the  first  system);  Chen  NNSC  Is  done 
once  for  each  additional  right-hand  side. 


PATH  - path  specification;  values  and  their  meanings  are  — 

1 perform  NROC,  NSFC,  and  NNFC. 

2 perform  NNFC  only  (NSFC  Is  assumed  to  hr.ve  been 

done  In  a manner  compatible  with  the  storage 
allocation  used  in  the  driver) . 

3 perform  NNSC  only  (NSFC  snd  NNFC  are  assumed  to 

have  been  done  In  a manner  compatible  with  the 
storage  allocation  used  in  Che  driver). 


1- 

n 


c 

1 

0 

No  Errors  Detected 

c 

1 

N4K 

Null  Row  in  A — Row  - K 

c 

1 

2N+K 

Duplicate  Entry  in  A — Row 

c 

1 

3N-HC 

Insufficient  Storage  In  NSFC 

c 

1 

AN+1 

Insufficient  Storage  In  NNFC 

c 

1 

5N-MC 

Null  Pivot  — Row  - K 

c 

1 

6N-fK 

Insufficient  Storage  In  NSFC 

c 

1 

7N+1 

Insufficient  Storage  in  NNFC 

G 

1 

8N-rtC 

Zero  Pivot  — Row  - K 

C 

1 

lON+1 

Insufficient  Storage  In  CDRV 

c 

1 

UN+1 

Illegal  PATH  Specification 

Various  errors  are  detected  by  the  driver  and  the  Individual 
subroutines. 


FLAG  - error  flag;  values  and  their  meanings  are  — 


- K 

— Row 


— Row 


Working  storage  Is  needed  for  the  factored  form  of  the  matrix 
H plus  various  temporary  vectors.  The  arrays  ISP  and  RSP  should  be 
the  same;  Integer  storage  Is  allocated  from  the  beginning  of  ISP 
and  real  storage  from  the  end  of  RSP. 

NSP  - declared  dimension  of  ISP  and  RSP;  NSP  generally  must 
be  larger  than  8N-f2  + 2K  (where  K - (number  of 
nonzero  entries  in  M)). 

nvlra  | ISP  - integer  working  storage  divided  up  Icto  various  arrays 
needed  by  the  subroutines;  ISP  and  RSP  should  be 
the  same  array. 

Size  • NSP. 

fvlra  I RSP  - real  working  storage  divided 
needed  by  the  subroutines; 
the  same  array. 

Size  - NSP. 

ESP  - If  sufficient  storage  was  available  to  perform  the 
symbolic  factorization  (NSFC),  then  ESP  is  set  to 
the  amount  of  excess  storage  provided  (negative  If 
Insufficient  storage  was  available  to  perform  the 
numeric  factorization  (NNFC)). 


up  Into  various  arrays 
ISP  and  RSP  should  be 


INTEGER  R(l),  C(l),  ICU),  IA{1),  JA(1),  ISP(l),  ESP,  PATH, 
FLAG,  TMP,  D,  Q,  U,  RMN,  ADD,  UMAX 
REAL  A(l),  B(l),  Z(l),  RSP(l) 


IF(PATH.LE.O  .OR.  PATH.GT.3)  GO  TO  111 
C******  Initialize  and  divide  up  temporary  storage 


******************* 


FLAG 

0 

IL 

1 

IJL 

IL 

+ N + 

1 

lU 

IJL 

+ N 

IJU 

lU 

+ N + 

1 

IRL 

IJU 

+ N 

JRL 

IRL 

+ N 

JL 

JRL 

+ N 

IRA 

NSP 

+ 1 - 

N 

D 

IRA 

JRA 

D 

- N 

TMP 

JRA 

Q 

TMP 

-(N  + 

1) 

JRU 

Q 

- N 

IRU 

JRU 

- N 

IF(JL 

GE. 

IRU) 

GO 

IF(PATH  .GT.  1) 

GO 

TO 


110 

10 


c 

c****** 


c****** 

10 


15 


20 


30 


Reorder  A If  necessary,  call 
RMN  - IRU  - JL 
ADD  - RMN/2 
JU  - JL  + ADD 
JLMAX  - ADD 
JUMAX  - RMN  - ADD 
DO  5 Il-l.N 

IF(C(II)  .NE.  II)  GO  TO  6 
CONTINUE 
GO  TO  7 


NSFC  if  flag  is  set  ************* 


CALL  NROC  (N,  1C,  lA,  JA,  A, 
IF(FLAG  .NE.  0)  GO  TO  100 


ISP(IL),  RSP(Q),  ISP(IU),  FLAG) 


CALL  NSFC 

(N,  R,  1C,  1A,JA,  JLMAX,  ISP(1L),1SP(JL),ISP(1JL),  JUMAX, 
1SP(IU),ISP(JU),1SP(IJU),  RSP(Q),  RSP(IRA),  RSP(JRA),  Z, 
ISP(IRL),ISP(JRL),  RSP(IRU),RSP(JRU),  FLAG) 

1F(FLAG  .NE,  0)  GO  TO  100 

See  if  enough  space  remains,  move  JU  next  to  JL  ************** 
JLMAX  - ISP(IJL-rtJ-l) 


JUMAX  - ISP(IJU-HI-l) 
LMAX  - ISP(IL+N)  - 1 
UMAX  - ISP(IU+N)  - 1 
IF(PATH  .GT.  1)  GO  TO 


20 


NEED  - JLMAX  + JUMAX  LMAX  + UMAX 
RMN  - RMN  + 3*N  f 1 
ESP  - RMN  - NEEi. 

IF(NEED  .GT.  RMN)  GO  TO  110 
JUOLD  - JU  - 1 
JU  - JL  + JLMAX  - 1 
IF  (JUMAX. LE.O)  GO  TO  20 
DO  15  Il-l, JUMAX 

1SP(JU+1I)  - ISP(JU0LD+1I) 

Call  remaining  subroutines 
JU  - JL  + JLMAX 
L - JU  + JUMAX 
U - L + LMAX 


*********************************** 


IF (PATH  .EQ 
CALL  NNFC 

(N,  R,  C,  IC, 


3)  GO  TO  30 


U,JA,A,  LMAX,  1SP(1I  ),ICr(JL),ISP(I 
RSP(D),  UMAX,ISP(1U),ISP(JU),ISP(IJU),RSP(U),  Z, 
RSP(TMP),  ISP(IRL),ISP(JRL),  FLAG) 

1F(FLAG  .NE.  0)  GO  TO  100 
RETURN 


).RSP(L), 

, d. 


CALL  NNSC 

(N,  R,  C,  1SP(1L),1SP(JL),ISP(1JL),RSP(L),  RSP(D),  ISP(IU), 
ISP(JU),ISP(1JU),RSP(U),  Z,  B,  RSP(TMP)) 

RETURN 


C **  ERROR:  Error  Detected  in  NROC,  NSFC,  NNFC,  or  NNSC 
100  RETURN 

C **  ERROR:  Insufficient  Storage 

110  FLAG  - 10*N  + 1 
RETURN 

C **  ERROR:  Illegal  PATH  Specification 

111  FLAG  - 11*N  + 1 
■ RETURN 

END 


r 


T 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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YALE  SPARSE  HATRIX  PACKAGE  - NONSYMMETRIC  CODES 
SOLVING  THE  SYSTEM  OF  EQUATIONS  Mx  - b 

I.  SUBROUTINE  NAMES 

Subroutine  names  and  functions  are  — 

(1)  NROC  for  reordering; 

(2)  NSFC  for  symbolic  factorization; 

(3)  NNFC  for  numeric  factorization  and  solution; 

(4)  NNSC  for  solution. 

II.  CALLING  SEQUENCES 

The  coefficient  matrix  can  be  processed  by  an  ordering  routine 
(e.g.,  to  reduce  fillln  or  ensure  numerical  stability)  before  using 
the  remaining  subroutines.  If  no  reordering  Is  done,  then  set 
R(I)  “ C(l)  “ IC(I)  ” 1 for  1=1, ...,N.  If  an  ordering  subroutine 
Is  used,  then  NROC  should  be  used  to  reorder  the  coefficient  matrix 
The  calling  sequence  Is  — 

( (matrix  ordering)) 

(NROC  (matrix  reordering)) 

NSFC  (symbolic  factorization  to  determine  where  fillln  will 
occur  during  numeric  factorization) 

NNFC  (numeric  factorization  into  product  LDU  of  unit  lower 
triangular  matrix  L,  diagonal  matrix  D,  and  unit 
upper  triangular  matrix  U,  and  solution  of  linear 
system) 

NNSC  (solution  of  linear  system  for  additional  right-hand 
side  using  LDU  factorization  from  NNFC) 

(If  only  one  system  of  equations  is  to  be  solved,  then  the 
subroutine  TRK  should  be  used.) 

III.  STORAGE  OF  SPARSE  MATRICES 

The  nonzero  entries  of  the  coefficient  matrix  M are  stored 
row-by-row  In  the  array  A.  To  Identify  the  Individual  nonzero 
entries  In  each  row,  we  need  to  know  in  which  column  each  entry 
lies.  The  column  Indices  which  correspond  to  the  nonzero  entries 
of  H are  stored  In  the  array  JA;  l.e..  If  A(K)  - M(I,J),  then 
JA(K)  • J.  In  addition,  we  need  to  know  where  each  row  starts  and 
how  long  It  Is.  The  Index  positions  in  JA  and  A where  the  rows  of 
M begin  are  stored  in  the  array  lA;  l.e..  If  M(I,J)  Is  the  first 
(leftmost)  entry  In  the  I-th  row  and  A(K)  • M(I,J),  then 
IA(I)  - K.  Moreover,  the  Index  In  JA  and  A of  the  first  location 
following  the  last  element  In  the  last  row  Is  stored  In  lA(N-fl). 
Thus,  the  number  of  entries  In  the  I-th  row  Is  given  by 
IA(I>1)  - IA(I),  the  nonzero  entries  of  the  I-th  row  are  stored 
consecutively  In 

A(IA(I)),  A(IA(I)+l),  ...,  A(1A(I+1)-1), 

and  the  corresponding  column  indices  are  stored  consecutively  In 
JA(IA(I)),  JA(IA(I)+1),  JA(IA(I+1)-1). 

For  example,  the  5 by  S matrix 

( 1.  0.  2.  0.  0.) 

( 0.  3.  0.  0.  0.) 

M - ( 0.  4.  5.  6.  0.) 

( 0.  0,  0.  7.  0.) 

( 0.  0.  0.  8.  9.) 
would  be  stored  as 

1123456789 


U I 1 3 4 7 8 10 
JA|1  32234445 
All  7.  I i S.  A.  7.  R Q. 


The  strict  upper  (lower)  triangular  portion  of  the  matrix 
U (L)  la  stored  In  a similar  fashion  using  the  arrays  lU,  JU,  U 
(IL,  JL,  L)  except  Chat  an  additional  array  lJU  (IJL)  Is  used  to 
compress  storage  of  JU  (JL)  by  allowing  some  of  Che  column  (cow) 
Indices  to  used  for  more  chan  one  row  (column)  (n.b.,  L Is  stored 
by  columns).  IJU(K)  (IJL(K))  points  to  the  starting  location  In 
JU  (JL)  of  entries  for  the  Kch  row  (column).  Compression  In  JU 
(JL)  occurs  In  two  ways.  First,  If  a row  (column)  1 was  merged 
Into  the  current  row  (column)  K,  and  the  number  of  elements  merged 
In  from  (the  tall  portion  of)  row  (column)  1 la  the  same  as  Che 
final  length  of  row  (column)  K,  then  the  Kth  row  (column)  and  Che 
call  of  row  (column)  I are  Identical  and  IJU(IC)  (IJL(K))  may  point 
CO  the  start  of  Che  call.  Second,  If  some  tall  portion  of  the 
(K-i)sC  row  (column)  Is  Identical  to  the  head  of  the  Kth  row 
(column),  then  IJU(K)  (IJL(K))  may  point  to  the  start  of  that  tail 
portion.  For  example,  the  nonzero  structure  of  the  strict  upper 
triangular  part  of  Che  matrix 
d 0 X X X 
0 d 0 X X 
0 0 d X 0 
0 0 0 d X 
0 0 0 0 d 

would  be  represented  as 

I 1 2 3 4 5 6 

lU  I 1 4 6 7 8 8 
JU  I 3 4 3 4 
IJU  I 1 2 4 3 

The  diagonal  entries  of  L and  U are  assumed  to  be  equal  to  one  and 
are  not  stored.  The  array  D contains  the  reciprocals  of  the 
diagonal  entries  of  the  matrix  D. 


C IV.  ADDITIONAL  STORAGE  SAVINGS 

C In  NSFC,  R and  1C  can  be  the  same  array  In  the  calling 

C sequence  If  no  reordering  of  the  coefficient  matrix  has  been  done. 

C In  NNFC,  2 and  ROW  can  be  the  same  array.  R,  C and  IC  can  all 

C be  the  same  array  if  no  reordering  has  been  done.  If  only  the 
C rows  have  been  reordered,  then  C and  IC  can  be  the  same  array. 

C If  the  row  and  column  orderings  are  the  same,  then  R and  C can  be 

C the  same  array. 

C In  NNSC,  R and  C can  be  the  same  array  If  no  reordering  has 

C been  done  or  If  the  row  and  column  orderings  are  the  same.  Z and  B 

C can  be  the  same  array;  however,  then  B will  be  destroyed. 

C 

C V.  PARAMETERS 

C Following  Is  a list  of  parameters  to  the  programs.  Names  are 

C uniform  among  Che  various  subroutines.  Class  abbreviations  are  — 

C n - INTEGER  variable 

C f - REAL  variable 

C V - supplies  a VALUE  to  a subroutine 

C r - returns  a RESULT  from  a subroutine 

C 1 - used  INTERNALly  by  a subroutine 

C a - ARRAY 

C 

C Class  I Parameter 

C 1 

C fva  I A - nonzero  entries  of  the  coefficient  matrix  M,  stored 
C I by  rows. 

C I Size  • number  of  nonzero  entries  In  H. 


C fva 


c 

1 

Size  - 

N. 

c 

nva 

1 

c 

- ordering 

of  the  columns  of  H. 

c 

1 

Size  - 

N. 

c 

fvra 

1 

0 

- reciprocals  of  the  diagonal  entries  of 

the  matrix 

c 

1 

Size  - 

N. 

c 

nr 

1 

FLAG 

- error  flag:  values  and  tiioir  meanings 

are  — 

c 

1 

0 

No  Errors  Detected 

c 

1 

N+K 

Null  Row  In  A — Row  « 

K 

c 

1 

2N+K 

Duplicate  Entry  in  A — 

Row 

- K 

c 

1 

3N4K 

Insufficient  Storage  for 

JL 

--  Row  - K 

c 

1 

4N+1 

Insufficient  Storage  for 

L 

c 

1 

3N-M1 

Null  Pivot  — Row  ■ K 

c 

1 

6N-MI 

Insufficient  Storage  for 

JU 

- Row  - K 

c 

1 

7N+1 

Insufficient  Storage  for 

U 

c 

1 

8M+K 

Zero  Pivot  — Row  - K 

c 

nva 

1 

lA 

- pointers 

to  delimit  the  rows  of  A 

c 

1 

Size  - 

N+1. 

C nvra 

C 

C 

C nvra 

C 

C 

C nvra 
C 

C nvra 
C 

C nva 
C 

C nvra 
C 

C nv 

C 

C 

C nvra 
L 

c;  nv 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


fvra 


nv 

nva 

fvra 


fra 


- right-hand  side  b. 


IJL  - pointers  to  the  first  element  In  each  column  in  JL, 
used  to  compress  storage  In  JL. 

Size  - N. 

lJU  - pointers  to  the  first  element  In  each  row  In  JU,  used 
to  compress  storage  In  JU. 

Size  • N. 

IL  - pointers  to  delimit  the  columns  of  L. 

Size  « N+1, 

lU  - pointers  to  delimit  the  rows  of  U. 

Size  - N+1. 

JA  - column  numbers  corresponding  to  the  elements  of  A. 

Size  * size  of  A. 

JL  - row  numbers  corct-spond  Ing  to  the  elements  of  L. 

Size  - JLMAX. 

JLMAX  - declared  dimension  of  JL;  JLIIAX  must  be  larger  than 
the  number  of  nonzeros  in  the  strict  lower  triangle 
of  M plus  fillin  minus  compression. 

JU  - column  numbers  corresponding  to  the  elements  of  U. 

Size  - JUHAX. 

JUMAX  - declared  dimension  of  JU ; JUMAX  must  he  larger  than 
the  number  of  nonzeros  in  the  strict  upper  triangle 
of  M plus  fillin  minus  compression. 

L - nonzero  entries  in  the  strict  lower  triangular  portion 
of  the  matrix  L,  stored  by  columns. 

Size  - LMAX. 

LMAX  - declared  dimension  of  L;  LMAX  must  be  larger  than 

the  number  of  nonzeros  in  the  strict  lower  triangle 
of  H plus  fillin  (1L(N+1)-1  after  NSFC). 

N - number  of  varlables/equatlons. 

R - ordering  of  the  rows  of  M. 

Size  ■ N. 

U - nonzero  entries  in  the  strict  upper  triangular  portion 
of  the  matrix  U,  stored  by  rows. 

Size  - UMAX. 

UMAX  - declared  dimension  of  U;  UMAX  must  be  larger  than 

the  number  of  nonzeros  in  the  strict  upper  triangle 
of  M plus  fillin  (IU(N+1)-1  after  NSFC). 

Z - solution  X. 

Size  ” N. 


Z 


C***  Subroutine  NROC 

C***  Reorders  rows  of  A.  leaving  row  order  unchanged 
C 

SUBROUTINE  NROC  (N,  1C,  lA,  JA.  A.  JAR,  AR,  P,  FLAG) 
C 

C Input  parameters:  N,  IC,  lA,  JA,  A 

C Output  parameters:  JA,  A,  FLAG 


Parameters  used  Internally: 

I  P - at  the  Kth  step,  P Is  a linked  list  of  the  reordered 
{ column  indices  of  the  Kth  row  of  A;  F(N-t-l)  points 

I to  the  first  entry  in  the  list. 

I Size  - N+1. 

1 JAR  - at  the  Kth  step, JAR  contains  the  elements  of  the 
I reordered  column  indices  of  A. 

I Size  - N. 

1 AR  - at  the  Kth  step,  AR  contains  the  elements  of  the 
I reordered  row  of  A. 

I Size  - N. 

INTEGER  IC(1),  U(l),  JA(1),  JAR(l),  P(l),  FLAG 
REAL  A(l),  AR{1) 


Q AAA***  For  each  nonempty  row  ******************************* 
DO  5 K-1,N 
JMIN  - IA(K) 

JHAX  " IA<K+[ ) - 1 
IF(JMIN  .GT.  .!)1AX)  GO  TO  5 
P(N+1)  -N  + 1 

C ******  Insert  each  element  in  the  list  ********************* 

DO  3 J-JMIN,JMAX 
NEWJ  - IC(JA(J)) 

I - N + 1 

1 1F(P(1)  .GE.  NEWJ)  GO  TO  2 

I - P<I) 

GO  TO  I 

2 IF(P(I)  .EQ,  NEWJ)  GO  TO  102 

P(NEWJ)  - P(l) 

P(I)  - NEWJ 
JAR  (NEWJ)  - JA(J) 

AR(NEWJ)  - A<J) 

3 CONTINUE 

C ******  Replace  old  row  in  JA  and  A ************************* 

1 - N + 1 
DO  4 J-JMIN,JMAX 
I - P(I) 

JA(J)  - JAR(l) 

4 A(J)  - AR(I) 

5 CONTINUE 
FLAG  - 0 
RETURN 

C 

C **  ERROR:  Duplicate  entry  in  A 
102  FLAG  = N + K 
RETURN 
END 
C 


! 


c 

c 

C***  Subroutine  NSFC 

C***  Symbolic  LDU-factorizatlon  of  nonsyrametrlc  sparse  matrix 


(compressed  pointer  storage) 

SUBROUTINE  NSFC 

(N,  R,  IC,  lA.JA,  JLMAX,  IL,  JL.  ML,  JUMAX,  lU,  JU , 1 JU , Q,  IRA, 
JRA,  IRAC,  IRL,JRL,  IRU,JRU,  FLAG) 


Input  variables: 
Output  variables: 


N,  R,  IC,  lA,  JA,  JLMAX,  JUMAX, 
IL,  JL,  IJL,  lU,  JU,  lJU,  FLAG. 


nla 


nla 
n la 
n la 


n la 
nla 


n la 
nla 


Parameters  used  Internally: 

0 - Suppose  M'  is  the  result  of  reordering  M.  If 

processing  of  the  Ith  row  of  M'  (hence  the  Ith 
row  of  U)  is  being  done,  0(J)  is  initially 
nonzero  if  M'(I,J)  is  nonzero  (J.GE.l).  Since 
values  need  not  be  stored,  each  entry  points  to  the 
next  nonzero  and  QfN+l)  points  to  the  first.  N+1 
indicates  the  end  of  the  list.  For  example,  if  N”9 
and  the  5th  row  of  M'  is 
OxxOxOOxO 
then  Q will  initially  be 

aaaa8aal03  (a-  arbitrary). 

As  the  algorithm  proceeds,  other  elements  of  Q 
arc  inserted  In  the  list  because  of  fillln. 

()  is  used  in  an  analogous  manner  to  compute  the 
Ith  column  of  L. 

Size  “ N+l. 

IRA,  - vectors  used  to  find  the  columns  of  M.  At  the  Kth 
JRA,  step  of  the  factor  ization,  IRAC(K)  points  to  Che 

IRAC  head  of  a linked  list  in  JRA  of  row  indices  I 

such  that  1 .GE.  K and  M(I.K)  is  nonzero.  Zero 
indicates  the  end  of  the  list.  IRA(I)  (l.GE.K) 
points  to  the  smallest  J such  that  J .GE.  X and 
M(I,J)  Is  nonzero. 

Size  of  each  = N. 

IRL,  - vectors  used  to  find  the  rows  of  L.  At  the  Kth  step 
JRL  of  the  factorization,  JRL(K)  points  to  the  head 

of  a linked  list  in  JRL  ot  column  Indices  J 
such  J .LT.  K and  L(K,J1  is  nonzero.  Zero 
indicates  the  end  of  the  list.  IRL(J)  (J.LT.K) 
points  to  Che  smallest  I sue):  that  1 .GE.  K and 
L(l , J ) is  nonzero. 

Size  of  each  “ N. 

IRU.  - vectors  used  in  a manner  analogous  to  IRL  and  JRL 
JkU  to  find  the  columns  of  U. 

Size  o£  each  “ N. 


Internal  variab'es: 

JLPTR  - points  to  the  last  position  used  in  JL. 

JUFTR  - points  to  the  last  position  used  in  JU. 

JMIN.JMAX  - are  the  indices  in  A or  U of  the  first  and  last 
elements  to  be  examined  in  a given  row. 

For  example,  JMIN-IA(K),  IMAX-1A(K+1  )-l , 

INTEGER  CEND,  (^,  REND,  Rl'.  VJ 

INTEGER  1A(1),  JA{1),  IRAM),  JRA(  1 ) . II. (1),  JL  ( 1 ) , IJL(l) 
INTEGER  IU(1),  JU(I),  IJUd),  IRL(M,  ,iRI.(l),  IRU(l),  JRU  ( M 
INTEGER  R(l),  IC(l),  Q(l).  IRAC(l),  I'LAi. 


c ******  Initialize  pointers  **************************************** 

NPl  - N + 1 
JLMIN  - 1 
JLPTR  - 0 
IL(1)  - 1 
JUMIN  - 1 
JUPTR  - 0 
IU(1)  - 1 
DO  1 K-1,N 
IRAC(K)  - 0 
JRA(K)  = 0 
JRL(K)  - 0 

1 JRU(IC)  - 0 

C ******  Initialize  column  pointers  for  A *************************** 

DO  2 K-1,N 
RK  - R(K) 
lAK  - lA(RK) 

IF  (lAK  .GE.  IA(RK+1))  GO  TO  101 
JAIAK  - IC(JA(1AK)) 

IF  (JAUK  .GT.  K)  GO  TO  105 
JRA(K)  - IRACfJAUK) 

IRAC(JAIAK)  - K 

2 IRA(K)  - lAK 
C 

C ******  For  each  column  of  E and  row  of  U ************************** 
DO  41  K-1,N 
C 

C ******  Initialize  Q for  computing  Kth  column  of  L ***************** 
Q(NPl)  - NPl 
LUK  - -1 

C ******  by  filling  in  Kth  column  of  A ****************************** 

VJ  - IRAC(K) 

IF  (VJ  .EQ.  0)  GO  TO  5 

3 QM  - NPl 

4 M - QM 

QM  - Q(M) 

IF  (QM  .LT.  VJ)  GO  TO  4 
IF  (QM  .EQ.  VJ)  GO  TO  102 
LUK  - LUK  + 1 
Q(M)  - VJ 
q(VJ)  - QM 
VJ  - JRA(VJ) 

IF  (VJ  .NE.  0)  GO  TO  3 

C ******  Link  through  JRU  ******************************************* 

5 LASTID  - 0 

LASTI  - 0 
IJL(K)  - JLPTR 
I - K 

6 I » JRU(I) 

IF  (I  .EQ.  0)  GO  TO  10 
QM  - NPl 
JMIN  - IRL(I) 

JMAX  - IJL(l)  + IL(I+1)  - IL(I)  - 1 
LONG  - JMAX  - JMIN 
JTMP  - JL(JMIN) 

IF  (JTMP  .NE.  K)  LONG  - LONG  + 1 
IF  (JTMP  .EQ.  K)  R(I)  - -R(I) 

IF  (LASTID  .GE.  LONG)  GO  TO  7 
LASTI  - I 
USTID  - LONG 
IF  (LONG  .LE.  0)  GO  TO  b 


7 


Q ****** 


Q ****** 
****** 

10 


Q ****** 


Q ****** 

a 


Q ****** 

15 


C ****** 


And  merge  the  corresponding  columns  Into  the  Kth  column  **** 
DO  9 J-JMIN,JMAX 
VJ  - JL(J) 

M - QM 
QM  - Q(M) 

IF  (QM  .LT.  VJ)  CO  TO  8 
IF  (QH  .EQ.  VJ)  CO  TO  9 
LUK  » LUK  + I 
Q(M)  - VJ 
Q(VJ)  - QM 
- VJ 

CONT INUE 
CO  TO  6 

LASTl  is  the  longest  column  merged  into  the  Kth  ************ 
See  if  It  equals  the  entire  Kth  column  ********************* 
QM  - Q(NPl) 

IF  (QM  .NE.  K)  CO  TO  105 
IF  (LUK  .EQ.  0)  CO  TO  17 
IF  (LASTID  .NE.  LUK)  GO  TO  11 

If  so,  JL  can  be  compressed  ******************************** 

IRLL  - IRL(LASTI) 

IJL(K)  - IRLL  + 1 

IF  (JL(IRLL)  .NE.  K)  IJL(K)  - IJL(K)  - 1 
GO  TO  17 

If  not,  see  if  Kth  column  can  overlap  the  previous  one  ***** 
IF  (JLMIN  .GT.  JLPTR)  GO  TO  15 
QM  - Q(QM) 

DO  12  J-JLMIN,  JLPTR 

IF  (JL(J)  - QM)  12,  13,  15 
CONTINUE 
GO  TO  15 
IJL(K)  - J 
DO  14  I^, JLPTR 

IF  (JL(I)  ,NE.  QM)  CO  TO  15 
QM  - q(qM) 

IF  (QM  .GT,  N)  GO  TO  17 
CONTINUE 
JLPTR  - J - 1 

Move  column  Indices  from  Q to  JL,  update  vectors  *********** 
JLMIN  - JLPTR  + 1 
IJL(K)  " JLMIN 
IF  (LUK  .EQ.  0)  GO  TO  17 
JLPTR  - JLPTR  + LUK 
IF  (JLPTR  .GT.  JLMAX)  GO  TO  103 
QM  - Q(NPl) 

DO  16  J-JLMIN,  JLPTR 
QM  - q(qM) 

JL(J)  - QM 
IRL(K)  - IJL(K) 

IL(K+1)  - IL(K)  + LUK 

Initialize  Q for  computing  Kth  row  of  U ******************** 
Q(NPl)  - NPl 
LUK  - -1 


C ******  by  filling  in  Kth  row  of  reordered  A *********************** 
RK  - R(K) 

JMIN  - IRA(K) 

JMAX  - IA(RK+1)  - 1 
IF  (JMIN  ,GT.  JMAX)  GO  TO  20 
DO  19  J-JMIN,JMAX 
VJ  - IC(JA(J)) 
qM  - NPl 

18  M - QM 
QM  - Q(M) 

IF  (QM  .LT.  VJ)  GO  TO  18 
IF  (QM  .EQ.  VJ)  GO  TO  102 
LUK  - LUK  + 1 
Q(M)  - VJ 
Q(VJ)  - QM 

19  CONTINUE 

C ******  Link  through  JRL,  ****************************************** 

20  LASTID  » 0 
LASTI  - 0 
IJU(K)  “ JUPTR 

I - K 

II  - JRL(K) 

21  I - II 

IF  (I  .EQ.  0)  GO  TO  26 
II  - JRL(I) 

QM  - NPl 
JMIN  - IRU(I) 

JMAX  - IJU(I)  + IU(I+1)  - IU(I)  - 1 
LONG  - JMAX  - JMIN 
JTMP  - JU(JMIN) 

IF  (JTMP  .EQ.  K)  GO  TO  22 

C ******  Update  IRL  and  JRL,  ***************************************** 

LONG  =■  LONG  + 1 

CEND  - IJL(I)  + IL(1+1)  - IL(I) 

IRL(I)  - IRL(I)  + 1 
IF  (IRL(I)  .GE.  CEND)  GO  TO  22 
J - JL(IRL(I)) 

JRL(I)  - JRL(J) 

JRL(J)  - I 

22  IF  (LASTID  .GE.  LONG)  GO  TO  23 

LASTI  - I 
LASTID  > LONG 

23  IF  (LONG  .LE.  0)  GO  TO  21 

C ******  And  merge  the  corresponding  rows  into  the  Kth  row  ********** 
DO  25  J-JMIN,JMAX 
VJ  - JU(J) 

24  M - QM 
QM  - Q(M) 

IF  (QM  .LT.  VJ)  GO  TO  24 

IF  (QM  .EQ.  VJ)  GO  TO  25 

LUK  - LUK  + 1 
Q(M)  - VJ 
Q(VJ)  - QM 
QM  - VJ 

25  CONTINUE 
GO  TO  21 

C ******  Update  JRL(K)  and  IRL(K)  *********************************** 

26  IF  (IL(K+1)  .LE.  IL(K))  GO  TO  27 

J - JL(IRL(K)) 

JRL(K)  - JRL(J) 

JRL(J)  - K 


! 


1 


l 

r 


Q *****:k 

27 


****** 


Q ****** 

28 


29 

30 


31 

C ****** 

32 


33 

3A 


C 

^ ****** 

35 


LASTl  Is  the  longest  row  merged  into  the  Kth  *************** 

See  If  It  equals  the  entire  Kth  row  ************************ 

QM  - Q(NPl) 

IF  (QM  .NE.  K)  GO  TO  105  i 

IF  (LUK  .EQ.  0)  GO  TO  34  ‘ ; 

IF  (LAST ID  .NE.  LUK)  GO  TO  28  i 

If  so,  JU  can  be  compressed  ******************************** 

IRUL  - I RU (LAST I)  j 

IJU(K)  - IRUL  + 1 ’ 

IF  (JU(IRUL)  .NE.  K)  IJU(K)  - IJU(K)  - 1 ;i 

GO  TO  34 

If  not,  see  if  Kth  tow  can  overlap  the  previous  one  ********  ; 

IF  (JUMIN  .GT.  JUPTR)  GO  TO  32  ‘ 

QM  - Q(qM)  I 

DO  29  J-JUMIN, JUPTR  ? 

IF  (JU(J)  - QM)  29,  30,  32  ] 

CONTINUE  ^ 

GO  TO  32  i 

IJU(K)  - J ) 

DO  31  I-J, JUPTR  I 

IF  (JU(I)  ,NE.  QM)  GO  TO  32  I 

QM  - Q(QM) 

IF  (QM  .GT,  N)  GO  TO  34 
CONTINUE 
JUPTR  - J - 1 

Move  row  Indices  from  Q to  JU,  update  vectors  ************** 

JUMIN  - JUPTR  + 1 
IJU(K)  - JUMIN 
IF  (LUK  .EQ.  0)  GO  TO  34 
JUPTR  - JUPTR  + LUK 
IF  (JUPTR  .GT.  JUMAX)  GO  TO  106 
QM  - q(NPl) 

DO  33  J-JUMIN, JUPTR 
QM  - Q(QM) 

JU(J)  - QM 
IRU(K)  - IJU(K) 

IU(K+1)  - IU(K)  + LUK 

Update  IRU,  JRU  *************************************■.  c 

I - K 

II  - JRU (I) 

IF  (R(I)  .LT.  0)  GO  TO  36 
REND  - IJU(I)  + IU(I+1)  - IU(I) 

IF  (IRU(I)  .GE.  RENO)  GO  TO  37 
J - JUdRUd)) 

JRU(I)  - JRU(J) 

JRU(J)  - I 
GO  TO  37 
Rd)  - -Rd) 

I - II 

IF  (I  .EQ.  0)  GO  TO  38 
IRUd)  - IRUd)  + 1 
GO  TO  35 


36 

37 


c 

Q ******  Update  IRA,  JRA,  IRAC  ************************************** 

38  I - IRAC(K) 

IF  (1  .EQ.  0)  GO  TO  41 

39  II  - JRA(I) 

IRA(l)  - IRA(l)  + 1 

IF  (IRA(I)  .GE.  U(R(1)+1))  GO  TO  40 
IRAI  - IRA(I) 

JAIRAl  - IC<JA(IRA1)) 

IF  (JAIRAl  .GT.  I)  GO  TO  40 
JRA(I)  - 1RAC( JAIRAl) 

1RAC( JAIRAl)  - I 

40  I - 11 

IF  (1  .NE.  0)  GO  TO  39 

41  CONTINUE 
C 

IJL(N)  - JLPTR 
IJU(N)  - JUPTR 
FLAG  - 0 
RETURN 
C 

C **  ERROR:  Null  Row  in  A 

101  FLAG  - N + RK 
RETURN 

C **  ERROR:  Duplicate  entry  in  A 

102  FLAG  - 2*N  + RK 
RETURN 

C 4*  ERROR:  Insufficient  Storage  for  JL 

103  FUG  - 3*N  + K 
RETURN 

C **  ERROR:  Null  pivot 

105  FUG  - 5*N  + K 
RETURN 

C **  ERROR:  Insufficient  Storage  for  JU 

106  FUG  - 6*N  + K 
RETURN 

END 

C 

C 

c 

c***  Subroutine  NNFC 

C***  Numerical  LDU-factorization  of  sparse  nonsymmetric  matrix  and 
C solution  of  system  of  linear  equations  (compressed  pointer 

C storage) 

C 

SUBROUTINE  NNFC 

* (N,  R,  C,  IC,  lA.JA.A,  LMAX.IL,  JL,IJL,L,  D,  UMAX,  lU,  JU,  IJU, 

* U,  Z.  B,  ROW,  IMP,  IRL,JRL.  FUG) 

C 

C Input  variables:  N,  R,  C,  1C,  lA,  JA,  A,  B,  IL,  JL,  IJL, 

C LMAX,  lU,  JU,  IJU,  UMAX 

C Output  variables:  Z,  L,  D,  U,  FUG 


Parameters  used  internally: 

I  IRL,  - vectors  used  to  find  the  rows  of  L.  At  Che  Kch  step 
I JRL  of  Che  factorization,  JKL(K)  points  to  Che  head 

I of  a linked  list  in  JRL  of  column  indices  J 

I such  J .LT.  K and  L(K,J)  is  nonzero.  Zero 

I indicates  the  end  of  the  list.  IRL(J)  (J.LT.K) 

I points  to  the  smallest  1 such  that  I .CL.  K and 

I L(1,J)  is  nonzero. 

I Size  of  each  « N. 

I ROU  - holds  intermediate  values  in  calculation  of  U and  L. 

I Size  • N. 

I TMP  - holds  new  right-hand  side  b'  for  solution  of  Che 
1 equation  Ux  ” b' . 

1 Size  « N. 


Internal  variables: 

JMIN,  JHAX  - indices  of  the  first  and  last  positions  in  a row  to 
be  examined. 

SUM  - used  in  calculating  TMP. 

INTEGER  RK.UMAX 
REAL  LKI 

INTEGER  R(l),  C(l),  IC(I),  IA(  I ) . JA(1),  1L(1),  JL(1),  IJL(l) 
INTEGER  IU(1),  JU(1),  UU(l),  IRL(l),  JRL(l),  FLAG 
REAL  A(l).  L(l),  D(l),  U(l),  Z(l),  B(l),  ROW(l),  TMP(l) 

******  Initialize  pointers  and  test  storage  *********************** 
IF(IL(N+1)-1  .GT.  LMAX)  GO  TO  104 
1F(IU(N+1)-1  .GT.  UMAX)  CO  TO  107 
DO  1 K-l.N 

IRL(K)  - IL(K) 

JRL(K)  - 0 
1 CONTINUE 


******  for  each  row  *********************************************** 

DO  19  K-1,N 

******  Reverse  JRL  and  zero  ROW  where  Kch  row  of  L will  fill  in  *** 
ROW(K)  - 0 
II  - 0 

IF  (JRL(K)  .EQ.  0)  GO  TO  3 

I - JRL(K) 

2 12  - JRL(l) 

JRL(I)  - II 

II  - I 
ROW(I)  - 0 
I - 12 

IF  <1  .NE.  0)  GO  TO  2 

******  Set  ROW  to  zero  where  U will  fill  in  *********************** 

3 JMIN  - IJU(K) 

JMAX  - JMIN  + IU(K+1)  - IU(K)  - 1 
IF  (JMIN  .GT.  JMAX)  GO  TO  5 


DO  4 J*JMIN,JMAX 
ROW(JU(J))  - 0 
Place  Kth  row  of  A in  ROW 
RK  - R(K) 

JMIN  - lA(RK) 

JMAX  - IA(RK+1)  - 1 
DO  6 J^MIN.JMAX 

R0W(IC(JA(J)))  - A(J) 
CONTINUE 


********************************** 


7 

****** 


Initialize  SUM,  and  link  through  JRL  *********************** 
SUM  - B(RK) 

I - II 

IF  (I  .Eq.  0)  CO  TO  10 

Assign  the  Kth  row  of  L and  adjust  ROW,  SUM  **************** 
LKI  - -ROW(I) 

If  U is  not  required,  then  comment  out  the  following  line  ** 
L(1RL(1))  - -LKI 
SUM  - SUM  + LKI  * TMP(I) 

JMIN  - IU(I) 

JMAX  - 1U(1+1)  - I 
IF  (JMIN  .GT.  JMAX)  GO  TO  9 
MU  - IJU(l)  - JMIN 
DO  8 J -JMIN,  JMAX 

R0W(JU(MU+J))  - ROW(JU(MU+J))  + LKI  * U (J ) 

I - JRL(I) 

IF  (I  .NE.  0)  CO  TO  7 


Assign  Kth  row  of  U and  diagonal  D,  set  TMP(K) 
IF  (ROW(K)  .EQ.  0)  CO  TO  108 
DK  - 1 / R0W<K) 

D(K)  - DK 
TMP(K)  - SUM  * DK 
IF  (K  .EQ.  N)  GO  TO  19 
JMIN  - 1U(K) 

JMAX  • IU(K+1)  - 1 
IF  (JMIN  .GT.  JMAX)  CO  TO  12 
MU  - IJU(K)  - JMIN 
DO  11  J-JMIN,JMAX 

U(J)  - R0W(JU(MU+J))  * DK 
CONTINUE 


Update  IRL  and  JRL,  keeping  JRL  In  decreasing  order  ******** 

I - II 

IF  (I  .EQ.  0)  GO  TO  18 
IRL(l)  - IRL(I)  + 1 

II  - JRL(I) 

IF  (IRL(I)  .GE.  IL(I+1))  GO  TO  17 
IJLB  = IRL(l)  - IL(I)  + IJL(I) 

J = JL(IJLB) 

IF  (I  ,GT.  JRL(J))  GO  TO  16 
J - JRL(J) 

GO  TO  15 
JRL(I)  - JRL(J) 

JRL(J)  - I 
I - II 

IF  (I  .NE.  0)  GO  TO  14 
IF  (IRL(K)  .GE.  IL(K+1))  GO  TO  19 
J - JL(1JL(K)) 

JRL(K)  - JRL(J) 

JRL(J)  - K 
CONTINUE 


r ' ^ 


c 

C ******  Solve  Ux  ” IMP  by  back  substitution  ********************** 

K - N 

DO  22  1-1, N 
SUM  - TMP(K1 
JMIN  = IU(K) 

JMAX  - 1U(K+1)  - 1 
IF  (JMIN  .GT.  JMAX)  GO  TO  21 
MU  - IJU(K)  - JMIN 
DO  20  J -JMIN,  JMAX 

20  SUM  - SUM  - U(J)  * TMP(JU(MU+J)) 

21  TMP(K)  - SUM 
Z(C(K))  - SUM 

22  K - K-1 
FLAG  - 0 
RETURN 

C 

C **  ERROR:  Insufficient  Storage  for  L 
104  FLAG  - 4*N  + 1 
RETURN 

C **  ERROR:  Insufficient  Storage  for  U 

107  FLAG  - 7*N  + 1 
RETURN 

C **  ERROR:  Zero  Pivot 

108  FLAG  - 8*N  + K 
RETURN 

END 

C 

c 

c 

C***  Subroutine  NNSC 

C***  Numerical  solution  of  sparse  nonsymmetrlc  system  of  linear 
C equations  given  LOU-factorlzatlon  (compressed  pointer  storage) 


SUBROUTINE  NNSC 


* (N,  R,  C,  IL,  JL,  IJL,  L,  D, 

lU, 

JU,  IJU,  U,  Z, 

B,  TMP) 

1,/ 

c 

Input  variables:  N,  R,  C,  IL, 

JL, 

IJL,  L,  D,  lU, 

JU,  IJU 

c 

Output  variables:  Z 

C 

Parameters  used  Internally: 

C fla  I TMP  - temporary  vector  vdilch  gets  result  of  solving  Ly  - b. 
C I Size  - N. 

C 

C Internal  variables: 

C JMIN,  JMAX  - Indices  of  the  first  and  last  positions  In  a row  of 

C U or  L to  be  used. 

C 

INTEGER  R{1),  C(l),  1L{1),  JL(l),  IJL(l),  1U(1),  JU  (I ) , IJU(l) 
REAL  L(l),  D(l),  U(l),  B(l),  Z(l),  TMP(l) 


c 

******  Set  IMP  to  reordered  B ************************************* 

DO  1 K-1,N 

1 TMP(K)  - B(R(K)) 

C ******  Solve  Ly  * b by  forward  substiCutLon  ********************** 

DO  3 K-l,N 
JMIN  - IL(K) 

JMAX  - IL(K+1)  - I 
TMPK  - -D(K)  * TMP(K) 

TMP(K)  - -TMPK 
IF  (JMIN  .(ir.  JMAX)  GO  TO  3 
ML  - IJL(K)  - JMIN 
DO  2 J-JMIN,JMAX 

2 TMP(JL(ML+J))  - TMP(JL(ML+J))  + TMPK  * L(J) 

3 CONTINUE 

C ******  Solve  Ux  • y by  back  subscicuclon  ************************* 

K - N 

DO  6 1-1, N 
SUM  - -TMP(K) 

JMIN  - 1U(K) 

JMAX  - 1U<K+1)  - 1 
IF  (JMIN  .GT.  JMAX)  GO  TO  5 
MU  - IJU(K)  - JMIN 
DO  4 J -JMIN,  JMAX 

4 SUM  - SUM  + U(J)  * TMP(JU(MU+J)) 

5 TMP(K)  - -SUM 
Z(C(K))  - -SUM 
K - K - I 

6 CONTINUE 
RETURN 
END 


Appendix  4 

Test  Driver  for  Sparse  Nonsymmetrlc  Matrix  Package 


7/31/77 


C 

C 

C 

C 

C 


C***  Program  NTST 

C***  Test  Driver  for  Nonsymmetrlc  Codes  in  Yale  Sparse  Matrix  Package 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


Variables : 

NG  - size  of  grid  used  to  generate  test  problem. 

N - number  of  variables  and  equations  (•  NC  x NG). 

lA  - INTEGER  onc-dimenslonal  array  used  to  store  row  pointers 

to  JA  and  A;  DIMENSION  - N+1. 

JA  - INTEGER  one-dimensional  array  used  to  store  column 

indices  of  nonzero  elements  of  M;  DIMENSION  “ number  of 
nonzero  entries  In  M. 

A - REAL  one-dimensional  array  used  to  store  nonzero  elements 
of  M;  DIMENSION  * ntsnber  of  nonzero  entries  in  M. 

X - REAL  one-dimensional  array  used  to  store  solution  x; 
DLMENStOS  - N. 

B - REAL  one-dlraenslonal  array  used  to  store  right-hand-side  b 
DLMENSION  - N. 

P - INTEGER  one- dimensional  array  used  to  store  permutation  of 
rows  and  columns  for  reordering  linear  system; 

DIMENSION  » N. 

IP  - INTEGER  one-dimensional  array  used  to  store  inverse  of 
permutation  stored  in  P;  DIMENSION  - N. 

NSP  - declared  dimension  of  one-dimenslonal  arrays  ISP  and  RSP. 

ISP  - INTEGER  one-dimensional  array  used  as  working  storage 
( equlvalenced  to  RSP);  DIMENSION  - NSP. 

RSP  - REAL  one- dimensional  array  used  as  working  storage 
(equlvalenced  to  ISP);  DIMENSION  - NSP. 

ESP  - INTEGER  amount  of  excess  storage  available 


INTEGER  lA(lOl).  JA(500).  P(IOO),  IP(IOO),  ISPdSOO),  ESP, 

* CASE,  PATH,  FLAG,  APTR,VP,VQ,  X,  XMIN,  UMAX,  Y,YMIN,YMAX 
REAL  A(500),  Z(IOO),  B(IOO).  RSP(I500),  NAME(3) 

EQUIVALENCE  (ISP(l),  RSP(l)) 

DATA  NSP/1500/,  EPS/lE-5/, 

* NAME(l)/'N'/.  NAME(2)/'T'/,  NAME(3)/'C'/ 

INDEX(I.J)  - NCn  + J - NG 


C 


NG  - 3 
N - NG*NG 


1 


f 


c 

c 

c 


******  CASE-1  ->  NDRV.  CASE-2  ->  TDRV,  CASE-3  ->  CDRV  *********** 
DO  5 CASE-1,3 

******  Set  up  matrix  for  five-point  finite  difference  operator  ***** 
APTR  - 1 
DO  2 I-1,NG 
DO  2 J-l.NG 

VP  - INDEX  (I,  J) 

P(VP)  - VP 
IP(VP)  - VP 
U(VP)  - APTR 
SUM  - 0 

XMIN  - MAXO  < 1.  I~l) 

XMAX  - MINO  (NC,  I+l ) 

VMIN  - MAXO  ( 1,  J~l) 

YMAX  - MINO  (NG,  J+1 ) 

DO  1 X-XMIN.XMAX 
DO  1 Y-YMIN,VMAX 

IF  ((X-I)  * (Y-J)  .NE.  0)  GO  TO  1 
VQ  - INDEX(X.  Y) 

JA(APTR)  - VQ 
A(APTR)  - 8 

IF  (VP  .LT.  VQ)  A(APTR)  - -1 
IF  (VP  .GT.  VQ)  A(APTR)  - -2 
SUM  - SUM  + A(APTR)  * VQ 
APTR  - APTR  + I 

1 CONTINUE 
B(VP)  - SUM 

2 CONTINUE 
IA(N+1  ) = APTR 
NZA  = IA(N+l)  - 1 


i 


C ******  Output  original  array  A ************************************ 

IF  (CASE.EQ.l)  PRINT  1001,  NG,NG 

1001  FORMAT  (/'  ***  FIVE-POINT  OPERATOR  ON  ', 

* II,  ' BY  ' II,  ' GRID  ') 

IF  (CASE.EQ.l)  PRINT  1002,  ( IA( I ) ,1-1 .N ) . LA(N+1) 

1002  FORMAT  (/'  COEFFICIENT  MATRIX:  '/ 

* /'  lA  (INDICES  OF  FIRST  ELEMENTS  IN  ROWS) 

* /(10I5)) 

IF  (CASE.EQ.l)  PRINT  1003,  (I, JA{I ) ,A(I ) , I-1,NZA) 

1003  FORMAT  (/'  JA  A 

* /'I  COLUMN  INDICES  MATRIX 

* /(I3,  110,  F16.5)) 

IF  (CASE.EQ.l)  PRINT  1004,  (B(I),  I-1,N) 

1004  FORMAT  RIGHT  HAND  SIDE  B: 

* /(5F10.5)) 


C ******  Call  ODRV  ************************************************** 

FUG  - 0 
PATH  - 1 
CALL  ODRV 

* (N,  IA,JA,A,  P,IP,  NSP,RSP,  PATH,  FUG) 

IF  (FUG.NE.O)  GO  TO  101 


I 

) 


T 


c 

Q ******  Output  orderlnR  of  variables/equations  ********************* 
IF  (CASE. eg. 1)  PRINT  1005,  (I, P(1 ) , IP (I) , I-l,N) 

1005  FORMAT  (/'  ROW/COLUMN  ORDERING  FROM  ODRV:  '/ 

* /'  P IP  ' 

* /'I  ROW/COL  ORDERING  INVERSE  ORDERING  ' 

* /(IJ,  no,  120)) 

C 

C ******  Call  NDRV  / TDRV  / CDRV  ************************************ 

PATH  - 1 


IF 

(CASE. Eg. 1) 

CALL  NDRV 

* 

(N,  P,P,IP, 

lA.JA.A,  B, 

X, 

NSP,  ISP.RSP.ESP, 

PATH, 

FLAG) 

IF 

(CASE. Eg.  2) 

CALL  TDRV 

* 

(N.  P,IP, 

lA.JA.A,  B. 

z. 

NSP,  ISP.RSP.ESP, 

FLAG) 

IF 

(CASE. Eg. 3) 

CALL  CDRV 

* 

(N,  P.P.IP, 

U,JA,A,  B, 

Z. 

NSP,  ISP.RSP.ESP, 

PATH, 

FLAG) 

IF 

( FLAG.  Eg.  0) 

GO  TO  3 

PRINT  1006,  NAME (CASE),  FLAG 

1006  FORMAT  (/'  ERROR  IN  ',  Al,  'DRV:  FLAG  - 15) 

CO  TO  5 

C 

C ******  Calculate  error  ******************************************** 

3 SUM  - 0 

DO  4 I-1,N 

4 SUM  - SUM  + ((Z<I)-I)/I)**2 
RMS  - SQRT(SUM/N) 

C 

C ******  Output  solution  and  error  measure  ************************** 

PRINT  1007,  NAME(CASE),  (Z(I),  I-1,N) 

1007  FORMAT  (/'  SOLUTION  FROM  ',  Al,  'DRV:  ' 

* /(5F10.5)) 

C 

IF  (RMS.LE.EPS)  PRINT  1008,  RMS 

1008  FORMAT  (/'  SOLUTION  CORRECT:  RMS  ERROR  - ',  1PE8, 2) 

IF  (RMS.GT.EPS)  PRINT  1009,  RMS 

1009  FORMAT  (/'  SOLUTION  INCORRECT:  RMS  ERROR  - ',  1PE8. 2) 

C 

PRINT  1010,  ESP 

1010  FORMAT  (/'  EXTRA  STORAGE  AVAILABLE  - ',  14) 

C 

5 CONTINUE 
STOP 

C 

C ******  Error  messages  ********************************************* 

101  PRINT  1013,  FLAG 

1013  FORMAT  (/'  ERROR  IN  ODRV:  FLAG  - ',  15) 

STOP 

END 


i 

i 


1 


1 


rl  II  HIM  I 

1 


Appendix  5 

Sample  .Qu tJ>Ht  From  Test  Driver 


***  FIVE-POINT  OPERATOR  ON  3 BY  3 GRID 
COEFFICIENT  MATRIX: 

lA  (INDICES  OF  FIRST  ELEMENTS  IN  ROWS) 


1 

4 

8 11  15  20 

JA  A 

I 

COLUMN 

INDICES  MATRIX 

1 

1 

8.00000 

2 

2 

-1.00000 

3 

4 

-1.00000 

4 

1 

-2.00000 

5 

2 

8.00000 

6 

3 

-1.00000 

7 

5 

-1.00000 

8 

2 

-2.00000 

9 

3 

8. 00000 

10 

6 

-1.00000 

U 

1 

-2.00000 

12 

4 

8.00000 

13 

5 

-1.00000 

14 

7 

-1.00000 

15 

2 

-2.00000 

16 

4 

-2.00000 

17 

5 

8.00000 

18 

6 

-1.00000 

19 

8 

-1. 00000 

20 

3 

-2.00000 

21 

5 

-2,00000 

22 

6 

8. 00000 

23 

9 

-1.00000 

24 

4 

-2.00000 

25 

7 

8. 00000 

26 

8 

-1 . 00000 

27 

5 

-2,00000 

28 

7 

-2.00000 

29 

8 

8.00000 

30 

9 

-1 . 00000 

31 

6 

-2.00000 

32 

8 

-2.00000 

33 

9 

8. 00000 

I 


1 


RIGHT  HAND  SIDE  B: 

2.00000  6.00000  14.00000  18.00000  14.00000 

23.00000  40.00000  31.00000  44.00000 

ROW/COLUMN  ORDERING  FROM  ODRV: 


P 

ROU/COL  ORDERING 
1 

3 

7 
9 
6 
5 
2 

4 

8 

SOLUTION  FROM  NDRV: 

1.00000  2.00000  3.00000 

6.00000  7.00000  8.00000 


IP 

INVERSE  ORDERING 
1 

7 
2 

8 
6 
5 

3 
9 

4 


4.00000  5.00000 

9.00000 


1 
1 

2 

3 

4 

5 

6 

7 

8 
9 


SOLUTION  CORRECT:  RMS  ERROR  " 6.36E-09 
EXTRA  STORAGE  AVAILABLE  = 1384 


SOLUTION  FROM  TDRV: 

1. 00000  2.00000  3.00000  4.00000  5.00000 

6.00000  7.00000  8.00000  9.00000 

SOLUTION  CORRECT:  RMS  ERROR  = 6.36E-09 
EXTRA  STORAGE  AVAILABLE  - 1412 
SOLUTION  FROM  CDRV: 

1.00000  2.00000  3.00000  4.00000  5.00000 

6.00000  7.00000  8.00000  9.00000 

SOLUTION  CORRECT:  RMS  ERROR  “ 6. 36E-09 


EXTRA  STORAGE  AVAILABLE  • 1364 


I 


